Final  Report  to 


European  Office  of  Aerospace  Research  and  Development 

(EOARD) 

Air  Force  Office  of  Scientific  Research 
Air  Force  Research  Laboratory 


Project  Title 

Development  of  a  robust  smart  antenna  array  signal-processing  algorithm 
based  on  antenna  array  auto-calibration  using  GPS  signals. 

Contract  No. 

F73001  F30602  99MV072 


Document  source 

Antenna  Research  Group,  ECE,  UL, 
Limerick  /  Dr.  M.  O’Droma 

Filing 

Date 

November  2002 

Document  ID  Version  No. 

ARG-ARFL-PJ  1.0/  VER  1 .7 

No.  of  pages 

74 

Antenna  Research  Group,  ECE,  UL 

ARFL 

Approved 

Prepared 

M. O’Droma 

J.  Mao 

A.  Goacher 

Report  Documentation  Page 

Form  Approved 

OMB  No.  0704-0188 

Public  reporting  burden  for  the  collection  of  information  is  estimated  to  average  1  hour  per  response,  including  the  time  for  reviewing  instructions,  searching  existing  data  sources,  gathering  and 
maintaining  the  data  needed,  and  completing  and  reviewing  the  collection  of  information.  Send  comments  regarding  this  burden  estimate  or  any  other  aspect  of  this  collection  of  information, 
including  suggestions  for  reducing  this  burden,  to  Washington  Headquarters  Services,  Directorate  for  Information  Operations  and  Reports,  1215  Jefferson  Davis  Highway,  Suite  1204,  Arlington 

VA  22202-4302.  Respondents  should  be  aware  that  notwithstanding  any  other  provision  of  law,  no  person  shall  be  subject  to  a  penalty  for  failing  to  comply  with  a  collection  of  information  if  it 
does  not  display  a  currently  valid  OMB  control  number. 

1 .  REPORT  DATE  2.  REPORT  TYPE 

01  OCT  2003  N/A 

3.  DATES  COVERED 

4.  TITLE  AND  SUBTITLE 

Development  of  a  robust  smart  antenna  array  signal-processing 
algorithmbased  on  antenna  array  auto-calibration  using  GPS  signals. 

5a.  CONTRACT  NUMBER 

5b.  GRANT  NUMBER 

5c.  PROGRAM  ELEMENT  NUMBER 

6.  AUTHOR(S) 

5d.  PROJECT  NUMBER 

5e.  TASK  NUMBER 

5f.  WORK  UNIT  NUMBER 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

University  of  Limerick  Limerick  Ireland 

8.  PERFORMING  ORGANIZATION 

REPORT  NUMBER 

9.  SPONSORING/MONITORING  AGENCY  NAME(S)  AND  ADDRESS(ES) 

10.  SPONSOR/MONITOR'S  ACRONYM(S) 

11.  SPONSOR/MONITOR'S  REPORT 
NUMBER(S) 

SPC  01-4075 

12.  DISTRIBUTION/AVAILABILITY  STATEMENT 

Approved  for  public  release,  distribution  unlimited 

13.  SUPPLEMENTARY  NOTES 

The  original  document  contains  color  images. 

14.  ABSTRACT 

15.  SUBJECT  TERMS 

16.  SECURITY  CLASSIFICATION  OF:  17.  LIMITATION  OF 

18.  NUMBER  19a.  NAME  OF 

a.  REPORT  b.  ABSTRACT  c.  THIS  PAGE  |J|J 

unclassified  unclassified  unclassified 

76 

Standard  Form  298  (Rev.  8-98) 

Prescribed  by  ANSI  Std  Z39-18 


Air  Force  Research  Laboratory  -EOARD  Contract  No.F73001  F30602  99MV072 


Air  Force  Research  Laboratory  -EOARD 


Contract  No.F73001  F30602  99MV072 


Table  of  Contents 

1 .  History . 5 

2.  Distribution  List  and  Acknowledgements . 5 

3 .  Referenced  Documents . 5 

4.  Glossary . 6 

5.  Introduction . 7 

5.1  General . 7 

5.2  Background . 7 

5.3  Description  of  Problem . 8 

5.4  Project  Objectives . 9 

6.  Global  Positioning  System  (GPS) . 10 

6.1  Introduction . 10 

6.2  GPS  Features . 1 1 

6.3  GPS  Signal  Characteristics . 12 

6.4  Using  GPS  to  Calibrate  Antenna  Arrays . 13 

7.  Antenna  Array  Model . 15 

7. 1  Modelling  the  Antenna  Array . 15 

7.2  Modelling  Systemic  Errors . 17 

8.  Parameter  Estimation  Methods . 20 

8.1  Introduction . 20 

8.2  Beamforming  Techniques . 20 

8.2.1  Coventional  Beamfonner . 21 

8.2.2  Constrained  Beamfonner  (Capon  Beamfonner) . 22 

8.3  Subspace  DOA  Estimation  Methods . 22 

8.3.1  Subspace  Properties . 22 

8.3.2  MUSIC  (Multiple  Signal  Classification) . 23 

8.3.3  ESPRIT  (Estimation  of  Signals  Parameters  via  Rotation  Invariance  Techniques)  .  23 

8.3.4  WSF  (Weighted  Subspace  Fitting) . 25 

8.4  Maximum  Likelihood  Estimation  Method  (MLEM) . 26 

8.5  Conclusion . 27 

9.  Sensitivity  Analyses  of  Parameter  Estimation  Methods . 28 

9. 1  Effects  of  Model  Enors  and  Sensitivity  Analysis  on  Subspace  Methods[26] . 28 

9.1.1  Effects  of  Model  Enors  on  the  Accuracy  of  DOA  Estimation . 29 

9. 1 .2  Effects  of  Model  Errors  on  Angular  Resolution . 30 

9. 1 .3  Conclusions  from  Sensitivity  Analyses . 31 

9.2  Effects  of  Model  Errors  and  Sensitivity  Analysis  on  Maximum  Likelihood  Methods  ..32 

9.2.1  The  Sensitivity  of  DOA  Estimate  Accuracy  to  Model  Errors . 33 

9.2.2  Failure  Threshold  of  ML  Algorithm . 34 

10.  Selection  of  Estimation  Method  and  Algorithm  Development . 36 

10.1  Selection  of  Estimation  Method . 36 

1 0.2  Algorithm  Development . 36 

10.2. 1  Maximum  Likelihood  Method . 36 

10.2.2  Iterative  MUSIC  Method . 40 

10.2.3  Separable  Sub-space  Method . 41 

U.  L.  Antenna  Research  Group  |  Page  3  of  75  |  Document  ID:  ARG-ARFL-PJ1.0/  VER  1.7 


Air  Force  Research  Laboratory  -EOARD 


Contract  No.F73001  F30602  99MV072 


10.2.4  Alternating  Iterative  Method . 45 

10.2.5  Global  Optimisation  Based  on  Simulated  Annealing . 48 

1 1 .  Simulation  Results . 52 

11.1  Convergence  Properties . 52 

1 1 .2  DO  A  and  Mutual  Coupling  Estimation  with  Unknown  Frequencies . 53 

11.3  Array  Sensor  Gain/Phase  Error  Calibration . 57 

11.3.1  Alternating  Iterative  Method . 57 

1 1 .3.2  Simulated  Annealing  Method . 58 

12.  Conclusions . 61 

12.1  Introduction . 61 

12.2  Antenna  array  model  and  systemic  errors . 61 

12.3  Parameter  estimation  methods . 62 

12.4  The  use  of  GPS  signals  to  estimate  the  MCM  and  gain/phase  errors . 62 

Appendix  A . 69 


U.  L.  Antenna  Research  Group 


Page  4  of  75 


Document  ID:  ARG-ARFL-PJ1.0/  VER  1.7 


Air  Force  Research  Laboratory  -EOARD 


Contract  No.F73001  F30602  99MV072 


1.  History 

Modified  by  Date  Version 

M.  O’Droma  19  Nov.2002  1.6 

M.  O’Droma/J.Mao  25  Nov.  2002  1.7 


Comments 

Edits 

Revision  of  some 

bibliographic 

infonnation 


2.  Distribution  List  and  Acknowledgements 

2.1  Distribution  List 

UL 

Dr.  M.  O  Droma 
Dr.  J.  Mao 
A.  Goacher 


EOARD 

Dr.  C.  Reuter 

AFRL/IFGC 

Dr.  K.  Kwiat 
R.  N.  Smith 


2.2  Acknowledgements: 

University  of  Limerick  would  like  to  acknowledge  the  fulsome  support  for  this  project  of 
European  Office  of  Aerospace  Research  and  Development  (EOARD),  Air  Force  Office  of 
Scientific  Research,  Air  Force  Research  Laboratory,  and  in  particular  the  individual  support  of 
Drs  Dan  McAuliffe,  John  C  Cleary,  Barry  McKinney,  Kevin  Kwiat,  Chris  Reuter  and  R.N. 
Smith. 

3.  Referenced  Documents 


Ref. 

No. 

Referenced  Document  Title 

Document  Source 

1 

Contract  No.  F73001  F30602  99MV072 

U.  L.  Antenna  Research  Group 


Page  5  of  75 


Document  ID:  ARG-ARFL-PJ1.0/  VER  1.7 


Air  Force  Research  Laboratory  -EOARD 


Contract  No.F73001  F30602  99MV072 


4.  Glossary. 


BER 

Bit  Error  Rate 

BPSK 

Binary  Phase-Shift  Keying 

CDMA 

Code  Division  Multiple  Access 

CRLB 

Cramer-Rao  Lower  Bound 

DGPS 

Differential  mode  GPS 

DOA 

Direction  of  Arrival 

DSP 

Digital  Signal  Processor 

ESPRIT 

Estimation  of  Signal  Parameters  via  Rotational  Invariance  Techniques 

EVD 

Eigenvalue  Decomposition 

GPS 

Global  Positioning  System 

GS 

Gram-Schmidt  (Orthoganalisation) 

MCM 

Mutual  Coupling  Matrix 

ML 

Maximum  Likelihood 

MSE 

Mean  Square  Error 

MUSIC 

Multiple  Signal  Classification 

MVDR 

Minimum  Variance  Distortionless  Response 

PRN 

Pseudo-Random  Noise 

RMSE 

Root  Mean  Square  Error 

SNR 

Signal-to-Noise  Ratio 

UAV 

Uninhabited  Air  Vehicle 

UCA 

Uniform  Circular  Array 

ULA 

Uniform  Linear  Array 
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5.  Introduction. 

5.1  General. 

For  the  U.S.  Air  Force  better,  more  rapid  decision-making  is  essential  to  achieving 
virtually  all  joint  war-fighting  capabilities  in  the  21st  century  battle-space.  War-fighters  must 
be  assured  the  capability  to  seamlessly  acquire,  store,  distribute  and  protect  their  information. 
Seamless  communications  span  the  globe;  interconnecting  command  echelons,  services,  and 
allies  worldwide  by  implementing  common  transport  protocols  and  dynamic  network 
management.  By  focusing  on  wide  bandwidth  capabilities  linked  to  the  Air  Force’s  current 
narrowband  tactical  systems,  including  mitigated  modems  to  recover  messages  during  nuclear 
and  naturally  disturbed  environments,  the  Air  Force  can  provide  the  correct  critical  information 
to  the  warrior  anywhere  in  the  world.  To  achieve  wideband  seamless  communications  it  is 
essential  to  use  smart  confonnal  wideband  (2  MFIz  -  2  GHz)  antennas.  This  has  the  added 
benefit  of  resulting  in  fewer  antennas  being  required  on  airplanes  and  Uninhabited  Air  Vehicles 
(UAVs). 

The  University  of  Limerick  has  ongoing  research  in  the  general  area  of  array  signal 
processing  algorithms  for  smart  antennas.  The  principal  general  goal  of  this  project  is  the 
development  of  improved  robust,  stable  and  accurate  algorithms  for  signal  source  bearing 
estimations  using  sampled  data  taken  from  the  sensors  of  phase  array  antennas  on  airplanes  and 
UAVs.  The  particular  goal  is  to  investigate  the  feasibility  of  using  Global  Positioning  System 
(GPS)  signals  to  facilitate  "on-the-fly"  auto-calibration  of  phase  array  antennas  and  to  develop 
a  means  of  doing  so. 


5.2  Background. 

While  the  potential  benefits  of  phase  array  antennas  are  well  known,  with  some  of  these 
having  been  exploited  in  various  applications  for  over  50  years,  the  developments  in  DSP 
hardware  and  software  are  now  making  it  feasible  to  exploit  more  and  more  of  these  benefits. 
This  has  led  to  a  renewed  interest  in  the  field  of  phase  array  antennas  and  to  an  upsurge  in  the 
development  of  parameter  estimation  algorithms.  Good  bibliographies  can  be  found  in  Krim 
&  Viberg's  review  paper  in  1996  [1]  and  in  Dr.  Tsoulos's  paper  [2], 


A  common  thread  running  through  the  published  literature  dealing  with  DSP  algorithm 
development  (some  of  these  algorithms  are  quite  ingenious,  innovative  and  unique)  is  the  use 


of  ideal  array  data  models.  This  avoids  consideration  =  the  more  practical  problems  of  real 


dynamic  array  system  characteristics  and  the  resultant  performance  degradation.  These 
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perfonnance-degrading  characteristics  are  mainly  the  sensor  gain/phase  uncertainty,  the  sensor 
location  uncertainty,  and  the  mutual  coupling  effect  between  the  sensor  array  elements.  In  this 
project  we  are  proposing  to  investigate  the  possibility  of  accounting  for  the  real  dynamic 
antenna  array  characteristics  by  means  of  an  "on-the-fly"  auto-calibration  system  using  GPS 
signals  to  produce  an  accurate  data  array  model.  If  this  is  achievable,  and  there  are  good 
reasons  to  believe  it  is,  then  signal  parameter  estimation  by  smart  antennas  could  become  a 
two-stage  process: 

(i)  auto-calibration  to  yield  an  accurate  array  data  model  and 

(ii)  the  signal  parameter  estimation  itself,  using  this  accurate  array  data  model,  with  one 
of  the  more  promising  modern  estimation  algorithms  (e.g.  maximum  likelihood  or 
subspace  methods). 

There  are  important  benefits  in  military  applications  for  this  technology  including  range 
extension,  capacity  enhancement,  higher  data  rates,  efficient  management  of  energy  resources 
and  better  BER  performance.  These  will  provide  greater  security  of  radio  communication 
channels  and  more  accurate  target  or  radio  source  pinpointing  and  tracking,  with  more  accurate 
control  of  antenna  pattern  nulls  and  peaks  -  an  important  attribute  in  jamming  counter¬ 
measures,  both  effecting  and  resisting  jamming. 

With  its  inherent  adaptive  capacity  to  enable  dynamic  narrow-cast  radio 
communications,  this  technology  can  contribute  towards  improvements  in  security,  immunity 
to  interference  and  robustness  of  communication  channels  whether  they  be  used  for 
telecommand,  telecontrol,  telemetry  or  multimedia  data/voice/video  exchange. 

5.3  Description  of  Problem. 

The  smartness  in  smart  antennas  lies  mainly  in  the  sophisticated  signal  processing 
algorithms.  These  have  the  task  of  extracting  the  required  signal  parameters  of  multiple  signal 
sources  from  the  data  received  at  the  antenna  array  sensors.  These  parameters  may  be  the 
number  of  signal  sources,  the  signal  frequencies,  the  directions  of  arrival,  polarisation,  etc.. 
This  is  a  non-trivial  multi-dimensional  estimation  problem,  and  thus  the  necessary  signal 
processing  algorithms  are  complex.  The  complexity  of  these  estimation  algorithms  is  further 
magnified  when  array  model  errors  such  as  the  inequality  of  sensor  channel  gain  and  phase 
characteristics,  sensor  location  errors,  and  the  presence  of  inter-sensor  mutual  coupling  are 
included  in  the  numerical  model. 
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One  way  of  reducing  the  complexity  of  the  estimation  algorithms  is  to  characterise  the 
array  model  more  accurately  by  estimating  the  net  effect  of  the  various  errors  and  using  this 
real  array  model  (as  opposed  to  the  ideal  array  model)  in  the  parameter  estimation.  This 
requires  signal  sources  with  known  parameters  (frequency  and  location)  covering  the  utilised 
array  manifold.  This  calibration  of  the  antenna  array  characteristics  will  lead  to  systems  that 
are  able  to  realise  more  reliable,  robust  and  accurate  parameter  estimation.  This  project 
proposes  to  take  an  important  step  in  this  direction  by  attempting  to  develop  a  capability  of  ‘on- 
the-fly’  auto-calibration  of  a  phase  array  antenna  system  using  the  known  characteristics  of 
GPS  signals  when  operating  in  an  otherwise  unknown  signal  environment. 

The  advantage  accruing  from  this  approach  is  that  whatever  the  parameters  are  which 
have  to  be  estimated  (directions  of  arrival,  signal  carrier  frequencies  etc.),  the  required  smart 
antenna  signal  processing  algorithms  can  now  use  an  accurate  array  data  model.  Thus  one 
significant  source  of  doubt  in  the  accuracy  of  the  estimated  parameters  is  removed.  At  the 
same  time,  removing  the  need  to  attempt  to  allow  for  unknown  changes  in  the  antenna  physical 
characteristics  since  the  ‘last’  calibration  reduces  the  complexity  of  the  estimation  algorithm. 

5.4  Project  Objectives. 

The  main  objective  of  this  project  was  to  develop  an  algorithm  that  will  use  the  known 
characteristics  of  GPS  signals,  which  are  received  by  the  phase  array  antenna,  to  execute  ’on- 
the-fly'  auto-calibration  of  the  antenna  array  thus  reducing  errors  due  to  unknown  sensor 
channel  characteristics,  array  characteristics  and  mutual  coupling  effects.  A  secondary 
objective  was  to  see  how  such  an  algorithm  could  then  be  combined  with  an  appropriate 
parameter  estimation  algorithm. 

A  successful  outcome  to  this  research  proposal  would  open  the  possibility  of  designing 
a  two-stage  smart  antenna  signal-processing  scenario  for  some  smart  antenna  applications: 

a)  "on-the-fly"  calibration  of  antenna  array  and  sensor  tuned  receiver  channels 

b)  estimation  of  the  required  unknown  receive  signal  parameters. 

The  result  should  be  that  the  accuracy  and  speed  of  parameter  estimation  would  be 
much  improved  over  that  of  existing  approaches. 
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6.  Global  Positioning  System  (GPS). 

6.1  Introduction. 

GPS  is  a  satellite  based  system  which  enables  the  positions  of  points  on  the  ground  or 
in  the  air  to  be  determined  with  high  accuracy  at  any  time  of  day  or  night  independent  of 
weather  conditions.  Over  the  past  years  GPS  has  gone  from  an  exciting  space-age  prospect  to 
an  invaluable  asset  in  many  fields  such  as  engineering  and  surveying  as  well  as  all  aspects  of 
navigation.  Currently  GPS  comprises  a  constellation  of  24  satellites.  There  are  4  satellites  in 
each  of  6  equispaced  orbits  that  have  an  inclination  of  55°  -  see  Fig  6.1. 


Fig  6. 1  Global  Positioning  System  (GPS)  Constellation. 

GPS  uses  the  principle  of  triangulation  to  calculate  a  user  receiver’s  position  relative  to 
a  set  of  known  points.  These  known  points  are  the  positions  of  the  satellites  (see  Fig  6.2)  as 
they  orbit  the  earth  and  broadcast  their  positions  and  a  consistent  time  standard  called  GPS 
System  Time.  The  receiver  uses  the  location  of  each  satellite,  the  system  time,  its  clock  biased 
relative  to  GPS  system  time  and  the  time  of  arrival  of  the  signal  to  compute  its  distance  from 
each  satellite.  More  sophisticated  user  equipment  also  uses  Doppler  readings  to  aid  in 
computing  the  receiver  velocity  or  to  provide  highly  accurate  navigation  solutions. 
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Fig  6.2  GPS  Constellation  -  planar  projection. 


6.2  GPS  Features. 

The  theory  behind  GPS  can  be  quite  easily  understood  without  going  into  the  complex 
mathematics  of  the  system.  The  distance  (range)  between  any  satellite  and  receiver  can  be 
calculated  by  measuring  the  time  taken  by  a  radio  signal  to  travel  from  the  satellite  to  the 
receiver.  The  geometry  of  the  system  is  such  that  if  three  ranges  are  measured  simultaneously 
then  the  user  can  calculate  his  position  using  a  computational  technique  known  as  resection. 
However,  in  order  to  solve  all  of  the  mathematical  unknowns  four  satellite  ranges  are  generally 
required.  These  satellites  must  be  visible  to  the  receiver  and  also  scattered  across  the  sky  in 
such  a  way  as  to  provide  a  good  geometric  arrangement  with  respect  to  the  receiver.  As  each 
satellite  transmits  mutually  exclusive  signals  it  is  possible  to  distinguish  between  information 
received  from  different  satellites.  Current  GPS  receivers  can  accept  information  from  up  to  8 
different  satellites  simultaneously.  The  type  of  infonnation  extracted  from  the  signal  and  the 
obtainable  accuracy  depends  on  the  type  of  receiver  used. 

There  are  two  categories  of  receiver  -  those  capable  of  reading  the  code  signal  only  and 
those  capable  of  reading  both  the  code  and  the  carrier  phase  signals  simultaneously.  Both 
types  of  receiver  support  two  unique  codes  -  the  C/A  (Coarse  Acquisition)  -code  and  the  P 
(Precise)  -code.  The  C/A-code  is  not  very  complex,  is  easy  for  a  receiver  to  lock  on  to  and  is 
available  to  all  civilian  users.  The  P-code,  on  the  other  hand,  although  more  difficult  to 
acquire,  is  more  accurate.  Code  only  receivers  can  be  used  in  one  of  two  modes  and  the 
accuracy  achieved  depends  on  the  method  of  use.  In  a  stand-alone  mode  only  one  receiver  is 
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necessary  in  order  to  compute  an  approximate  position  by  resection.  The  positional  accuracy 
is  generally  good  to  ±20m,  but  it  is  only  reliable  to  ±100m.  The  American  military,  who 
control  the  GPS  satellites,  have  the  ability  to  downgrade  the  satellite  signal  through  a  process 
known  as  Selective  Availability  (SA)  and  it  is  this  that  limits  the  accuracy  of  code  only 
receivers.  Although  SA  degrades  the  detennination  of  absolute  positioning,  it  has  little  effect 
on  the  relative  position  between  two  close  receivers  as  the  error  in  the  received  signals  will  be 
similar  in  both  receivers  and  thus  cancel  out.  Accuracies  of  ±lm  can  be  obtained  by  using 
two  or  more  code-only  receivers  operating  simultaneously;  this  is  known  as  differential  mode 
GPS  (DGPS).  In  this  mode  of  operation  the  resection  solutions  at  each  ground  station  can  be 
related  to  each  other  producing  accurate  co-ordinate  differences  between  the  two  stations  rather 
than  absolute  positions. 

Since  metre  level  accuracy  is  inadequate  for  most  surveying  or  engineering  jobs  a 
second  category  of  receivers,  known  as  code  and  carrier  phase  receivers,  must  be  used  to 
achieve  the  required  accuracy.  These  receivers  are  only  designed  for  use  in  differential  mode. 
They  record  both  the  C/A-code  and  the  P-code  (when  available)  and  also  use  the  carrier  phase 
as  a  means  of  fine  interpolation.  Each  satellite  transmits  two  carrier  frequencies,  LI  and  L2, 
and  consequently  there  are  two  types  of  code  and  phase  receiver  available  -  single  frequency, 
which  use  the  information  available  on  LI  only,  and  dual  frequency,  which  use  both 
frequencies. 

6.3  GPS  Signal  Characteristics. 

GPS  satellite  transmissions  utilise  direct  sequence,  spread  spectrum  (DSSS) 
modulation.  DSSS  provides  the  structure  for  the  transmission  of  ranging  signals  and  essential 
navigation  data  such  as  satellite  ephemeredes  and  satellite  height.  The  ranging  signals  are 
pseudo-random  noise  (PRN)  codes  that  binary  phase  shift  key  (BPSK)  modulate  the  satellite 
carrier  frequencies.  These  codes  look  like  and  have  spectral  properties  similar  to  random 
binary  sequences,  but  are  actually  detenninistic.  They  have  a  predictable  pattern  that  is 
periodic  and  can  be  replicated  by  a  suitably  equipped  receiver.  Each  GPS  satellite  broadcasts 
two  PRN  ranging  codes;  a  ‘short’  C/A-code  and  a  ‘long’  P-code.  The  C/A-code  has  a  period 
of  1  msec  and  repeats  constantly,  whereas  the  P-code  is  a  seven-day  sequence  that  repeats  every 
Saturday/Sunday  midnight.  At  the  time  of  writing  the  P-code  is  encrypted  and  is  known  as 
the  Y-code.  This  Y-code  is  only  accessible  to  selected  users  through  the  use  of  cryptography. 
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The  GPS  satellites  transmit  on  two  carrier  frequencies  called  LI,  the  primary 
frequency,  and  L2,  the  secondary  frequency.  Each  carrier  frequency  is  modulated  by  the 
navigation  data  message  and  then  the  spectrum  spread  by  modulation  with  a  PRN  sequence 
unique  to  the  satellite.  All  satellites  transmit  on  the  same  two  carrier  frequencies,  but  their 
signals  to  not  interfere  significantly  because  each  has  a  unique  PRN  code  and  the  PRN  codes 
are  selected  to  be,  as  far  as  possible,  uncorrelated  with  each  other.  Because  the  PRN  codes  are 
nearly  uncorrelated,  the  signals  from  each  satellite  can  be  detected  and  separated  using  code 
division  multiple  access  (CDMA).  In  order  to  track  one  satellite  in  common  view  with 
several  others  using  the  CDMA  technique,  the  receiver  must  replicate  the  PRN  sequence  for 
the  desired  satellite  along  with  the  replica  carrier  signal,  including  Doppler  effects. 

Two  carrier  frequencies  are  transmitted  to  enable  two-frequency  users  to  measure  the 
ionospheric  delay.  This  delay  effects  the  accuracy  of  the  time  measurement  and  can  be 
measured  using  two  frequencies  since  this  delay  is  related  by  a  scale  factor  to  the  difference  in 
signal  time-of-arrival  for  the  two  carrier  frequencies.  Single-frequency  (LI  only)  users  must 
estimate  the  ionospheric  delay  using  modeling  parameters  that  are  broadcast  to  the  user  in  the 
navigation  message. 

6.4  Using  GPS  to  Calibrate  Antenna  Arrays. 

GPS  is  a  system  that  enables  a  user  to  calculate  his  position  from  the  times-of-arrival  of 
the  signals  from  a  number  of  satellites.  It  is  a  simple  step,  knowing  the  location  of  the  user 
and  the  positions  of  the  various  satellites  in  space,  to  calculate  the  directions-of-arrival  (DOAs) 
of  the  satellite  signals.  It  is  the  knowledge  of  these  DOAs,  in  azimuth  and  elevation  with 
respect  to  the  reference  plane  of  the  antenna  array,  that  will  enable  the  calibration  of  the 
antenna  array  and  thus  reduce  errors  due  to  unknowns  such  as  physical  distortions  and  mutual 
coupling  effects.  Fig.  6.3  shows  how  four  GPS  satellites  could  be  used  as  calibration  sources 
with  known  azimuth  and  elevation. 
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Fig.  6.3  GPS  Signals  as  Calibration  Sources 
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7.  Antenna  Array  Model. 

7. 1  Modelling  the  Antenna  Array. 

Consider  an  antenna  array  composed  of  M  antenna  sensors  arbitrarily  located  in  space 
and  assume  that  a  signal  s(t)  impinges  on  the  array.  Let  xm  be  the  propagation  time  delay  of 
the  signal  at  the  111th  sensor,  related  to  a  fixed  point.  Assuming  that  the  antenna  can  be 
modelled  as  a  linear  system  and  letting  hm(t)  be  the  impulse  response  of  the  111th  sensor,  then  the 
output  of  the  111th  sensor  can  be  written  as: 

Xm(t)  —  hm(t)  oc  s(t  -  Tm)  +  nra(t) 

where  adenotes  convolution  and  nm(t)  is  the  additive  noise,  which  is  independent  of  the 

signal. 

Let  signal  s(t)  be  a  modulated  signal  with  carrier  angular  frequency  coc  (=  27tfc)  and 
modulation  amplitude  a(t)  and  phase  cf)(t)  such  that: 

s(t)  =  a(t)  eos(coct  +  0(t)) 

If  the  signal  being  considered  is  narrowband,  i.e.  the  envelope  variations,  amplitude  a(t) 
and  phase  c[)(t),  of  the  signal  vary  slowly  relative  to  the  propagation  time  across  the  array,  then: 

a(t  -  v)  =  a(t),  cf)(t  -  v)  =  0(0 

Hence,  we  may  write  the  received  signal  s(t)  as: 

s(t-  V)  —  a(t  -  Tm)  eos(coc(t  -  v)  "F  0(t  ■  xm)) 

=  a(t)  eos(coc(t  -  xm)  +  0(t)) 

and  the  output  of  the  mth  sensor  may  be  written  as: 

Xm(t)  hm(t)  oc  s(t  -  Xm)  T  nm(t) 

=  hm(t)  a  a(t)  eos(coc(t  -  xm)  +  0(t))  +  nm(t) 

The  narrowband  assumption  is  most  often  satisfied  when  the  signal  bandwidth  is  much 
smaller  than  the  carrier  frequency  and  the  signal  propagation  time  across  the  array  is  smaller 
than  the  inverse  of  the  array  aperture.  This  assumption  can  also  be  satisfied  for  wideband 
signals,  such  as  in  the  case  of  CDMA,  if  the  frequency  response  of  each  antenna  sensor  is 
approximately  flat  over  the  signal  bandwidth.  Using  this  narrowband  assumption,  the 
propagation  time  delay  across  the  array  can  be  modelled  as  a  simple  phase  shift  and  the  signals 
have  bandpass  characteristics  about  the  centre  frequency.  The  complex  representation  of  the 
output  of  the  mth  sensor  can  then  be  given  as: 

Xm(t)  Hm(coc)  exp(-jcocxm)  s(t)  nm(t) 
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where  Hm(coc)  is  the  frequency  response  of  the  mth  sensor  (i.e.  the  Fourier  transform  of 

hm(t)). 

Given  the  assumption  that  each  antenna  sensor  can  be  modelled  as  a  linear  element,  the 
resultant  sensor  output  signal  is  a  superposition  of  the  p  individual  signals  received.  Here,  we 
introduce  the  parameter  vector  r\„  (n  =  1,2,. .  ,,p)  to  denote  the  collection  of  parameters,  such  as 
bearing,  elevation,  frequency,  polarisation  and  so  on,  associated  with  the  nth  signal.  The 
impulse  response  and  time  delay  of  each  sensor  can  be  represented  as  a  function  of  the 
parameter  vector  r|n  so  that  the  data  model  for  the  mth  sensor  output  can  be  rewritten  as: 

Xm(t)=EHm(tli)eXP(_ja)cTm(Tl1))Si(t)+nm(t)=Eam(tli)Si(t)  +  nm(t) 
i=l  i=l 

Since  the  array  comprises  M  sensors,  the  output  vector  for  the  array  is  given  by: 

fXl(tn  fa^n  km 

x(t)=  :  =2  :  Sl(t)+  ; 

VXM(0y  a  M  (^li  ) y  vnM(0y 

=  Za(1l.)s.(t)  +  n(t) 

=  [a (il i  ),*•',  a (r) p )  :  +n(t) 

Ml 

=  A(ii)s(t)  +  n(t) 

where  the  vector  a(ri)  is  referred  to  as  the  array  response  vector  or  steering  vector  and 
is  given  by : 

a(Tl) =  [H, (n)exp(- j coc x j (r|))  H2(t])exp(- j(Ocx2 (i|))  HM(n)exp(- jcocxM(n))]T 

and  a(ri)  is  a  function  of  the  parameter  vector  r),  r)  =  [  r)  iT,  . . . ,  r)pT]T. 

If  there  are  N  (<M)  elements  (different  parameters)  in  a(ri),  then  r)  will  trace  an  N- 
dimensional  surface  in  CM  as  p  is  varied  over  the  parameter  space.  This  surface  is  referred  to 
as  the  array  manifold  and  is  denoted: 

A  =  {a(r|)  :  r|  e  0} 

where  0  denotes  the  set  of  all  possible  parameter  vectors  of  interest. 

The  noise  process,  n(t),  is  assumed  to  be  a  zero-mean,  stationary  random  process  and 
the  probability  of  distribution  of  the  noise  is  assumed  to  be  circular  complex  Gaussian,  i.e.  the 
real  and  imaginary  parts  of  n(t)  are  independent,  and  with  the  second  order  moments: 

E{n(t)nH(x)}  =  Q  8(t,x) 

E{n(t)nT(x)}  =  0 
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where  E{.}  and  8  represent  the  statistical  expectation  and  Kronecker  delta  function 
respectively.  In  this  analysis,  the  background  noise  is  assumed  to  be  spatially  white  so  that 
the  noise  covariance  matrix,  Q,  is  a  scaled  identity  matrix  a  I,  where  a  is  the  noise  variance. 


In  this  paper,  since  we  are  considering  signals  of  known  frequency,  we  are  only 
interested  in  the  parameters  of  signal  azimuth  and  elevation.  The  array  manifold  can, 


therefore,  be  described  as: 


A  =  {a(0,<])) :  0,4>  e  0} 


which  is  a  two-dimensional  surface  in  space  CM,  where  0  and  0  denote  the  signal 
azimuth  and  elevation  respectively. 

Thus,  the  output  vector  for  the  array  can  be  rewritten  as: 

x(t)  =  Xa(0i^i)si(t)+n(t) 

i=l 

=  A(0,<p)S(t)  +  n(t) 

where  the  steering  vector: 

a(®i.<l)i )=  [H](6i,(|>i)e_jc*tXl*0|'*i’  H2(9i,(|>1)e~jc*‘T;(0i4,)  -  *>]T 

If  we  assume  each  array  sensor  has  the  same  response,  then  Hi(0n,(f)n)  =  HM(0n,(|>n)  = 
H(0n,c|)n).  Thus,  a(0i,cf>i)  is  a  function  of  xm  only  and  can  be  written: 


a(0,,0,)  =  H(0i,^i)[e'j“cXl(e'4i)  e‘j“^(eiA) 


In  the  ideal  case,  the  time  delay  at  each  element  can  be  calculated  for  any  0  and  0  from 
the  array  geometry  and  the  signal  frequency,  and  the  calibration  of  the  array  depends  only  on 
detennining  the  sensor  response  H(0,(f>). 


7.2  Modelling  Systemic  Errors. 

Whilst,  in  the  ideal  case,  the  steering  vector  can  be  calculated  from  knowledge  of  the 
array  geometry,  in  practice  there  are  a  number  of  sources  of  error  in  any  system.  These  may 
be  grouped  under  two  headings. 

Firstly,  there  are  errors  in  amplitude  and  phase  in  the  tracking  between  the  different 
sensors  and  their  associated  circuitry.  These  may  be  due  to  inherent  imbalances  between  the 
sensors  and  their  associated  circuitry  (amplifiers,  cables,  digitisers,  etc.)  or  they  may  be 
changes  in  any  part  of  the  system  due  to  temperature  or  aging  effects. 
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Consider  p  radiating  sources  observed  by  an  array  of  M  antenna  sensors.  The  signal 
output  of  the  mth  sensor  can  be  described  by: 

p 

Xm  (t)  =  X  amS„  (t  -  Lnn  “  ¥,n  )  +  \  (0 
n=l 

where  sn(t)  (n  =  1,2,  . .  .p)  represents  the  received  signals,  0Cm  and  \|/m  are  the  gain  and 
phase  delay  associated  with  the  mth  sensor  and  its  circuitry  and  xmn  is  the  delay  relative  to  a 
reference  point  associated  with  the  signal  propagation  from  the  nth  source  to  the  mth  sensor. 

In  matrix  notation  the  array  output  can  be  expressed  as: 

x(t)  =  rA(e,<p)s(t)+N(t) 

where  - 

X(t)=[xi(t)  x2(t)  xm(0]T 
S(t)  =  [s, (t)  s2(t)  •••  sM(t)]T 
N(t)=[n1(t)  n2(t)  •••  nM  (t)]T 

r  =  diag[a1ej°Wl  a2eJ“c'1'2  •  •  •  aVIc  |0)c'|,M  ]'  contains  the  amplitude  and  phase  errors 

A(0,(p)=[a(0,,01)  a(e2,4>2)  •••  a(0N,(|)N )] 

Secondly,  there  are  mutual  coupling  effects  between  the  antenna  sensors  that  make  up 
the  array.  When  antenna  sensors  are  in  close  proximity,  typically  less  than  half  a  wavelength, 
the  impedance  and  polar  response  of  each  sensor  is  affected  by  the  electromagnetic  coupling  to 
its  adjacent  sensors.  This  will  distort  the  radiation  pattern  of  the  array  and  hence  modify  the 
array  manifold.  The  output  of  the  mth  sensor,  in  the  absence  of  other  sources  of  error,  may 
now  be  written  as: 

M 

X  m  (t )  =  X  C  a  ■  (e,  <t>)s(t) + N(t ) 

i=l 

where  Cmj  (i,  m  =  1,2,  . .  .M)  is  the  mutual  coupling  factor  and  models  the  mutual 
coupling  effect  of  the  ith  sensor  on  the  mth  sensor  and  ai(0,c[))  is  the  ith  row  of  the  steering  matrix 

a(0,4>). 

Thus,  the  output  of  the  array  is  given  by: 

X(t)  =  C  A(0,0)  S(t)  +  N(t) 

where  C  is  an  M  x  M  complex  matrix  taking  account  of  all  effects  due  to  mutual 
coupling  between  sensors  and  is  referred  to  as  the  Mutual  Coupling  Matrix.  In  general,  the 
matrix  C  has  no  special  structure;  however,  if  the  array  is  unifonn,  then  the  matrix  will  be 
structured  because  coupling  between  adjacent  sensors  is  common  to  sensors  occupying  a 
similar  position  in  the  array  and  coupling  between  non-adjacent  sensors  may  be  ignored. 
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When  both  types  of  error  are  considered  the  array  output  may  be  written  as: 

X(t)  =  T  C  A(0,<|>)  S(t)  +  N(t) 

=  A’(0,(|))  S(t)  +  N(t) 

where  A’  represents  the  actual  steering  matrix,  allowing  for  all  sources  of  error,  at 
frequency  fc.  It  is  this  matrix  that  will  be  calculated  from  GPS  measurements  and  provide  a 
calibration  of  the  antenna  array. 
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8.  Parameter  Estimation  Methods. 

8.1  Introduction. 

As  shown  in  section  7.1,  the  output  vector  for  a  perfect  array  (i.e.  ignoring  amplitude  or 
phase  imbalances  and  the  effects  of  mutual  coupling  between  array  sensors)  can  be  given  by: 

X(t)  =  A(6,d)  S(t)  +  N(t) 

where  the  array  manifold  A(0,c|))  can  be  calculated  from  the  array  geometry  and  the 
operating  frequency. 

This  means  that,  assuming  the  array  manifold  is  unambiguous,  it  is  possible,  by  suitable 
processing  of  the  received  signal,  to  calculate  the  azimuth  and  elevation  of  the  signal.  This 
involves  using  spatial  processing  technology  to  generate  a  spatial  data  matrix  from  which  the 
desired  parameters  can  be  estimated  by  using  an  appropriate  algorithm. 

In  this  section  we  will  review  some  of  the  main  estimation  methods;  beamforming, 
subspace-based  methods  and  maximum  likelihood  methods. 

8.2  Beamforming  Techniques. 

Digital  beamforming  is  a  product  of  merging  antenna  technology  and  digital  signal 
processing  technology.  Early  applications  of  digital  beamfonning  were  developed  for  sonar 
and  radar  systems.  One  of  their  major  advantages  is  that  all  of  the  desired  information  being 
carried  by  the  signals  may  be  extracted  in  the  fonn  of  digital  streams  and  made  available  for 
processing  in  the  beamfonner.  Using  this  infonnation  and  suitable  algorithms,  the 
beamfonner  can  be  driven  to  produce  different  types  of  beams  such  as  scanning  beams, 
multiple  beams,  shaped  beams  or  beams  with  steered  nulls.  There  are  many  different 
configurations  that  may  be  used  to  achieve  these  beamformins  and  Fig.  8.1  depicts  a  typical 
beamfonner  structure.  The  beamfonner  forms  a  linear  combination  of  the  antenna  outputs, 
first  multiplying  each  output  by  a  complex  weight  and  then  summing  them  together. 
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Fig.  8.1  A  typical  beamformer  structure. 


8.2.1  Coventional  Beamformer. 


The  conventional  beamformer  is  a  natural  extension  of  classical  Fourier-based 
spectral  analysis  [3]  to  array  data.  The  beamformer  maximizes  the  power  of  the  array  output 
for  a  given  input  signal: 

x(t)  =  a(e,4>)s(t)  +  n(t) 

If  we  want  to  maximize  the  output  power  from  a  specific  direction  (0,(|))  then,  for  the 
conventional  beamfonner,  the  output  of  the  beamfonner  at  time  t,  y(t),  is  given  by  a  linear 
combination  of  the  data  at  M  antenna  sensors  at  time  t: 


M 

y(t)  =  ^w»Xm(t) 

m=l 


where  xm(t)  is  the  signal  from  the  mth  antenna  sensor  of  the  array  and  wm  is  the  weight 
applied  to  xm(t).  In  vector  fonn,  this  can  be  written  as: 

y(t)  =  wHx(t) 

Hence,  the  problem  of  maximizing  the  output  power  is  formulated  as: 
argmaxP(w)  =  argmaxE{wHx(t)xH(t)w} 

w  w 

=  argmaxwHE{x(t)xH(t)}w 


=  arg  max 

w 


{E|s(t)|2  |wHa(0,  cj>) 


2  2  1  |2  -v 

+  0  w  } 


Here  we  assume  the  noise  is  spatially  white  and  the  norm  of  w  is  constrained  to  |w|=l. 
Then,  the  resulting  solution  is: 

a(0,4>) 


wBF  = 


VaH(0,0)a(0,0) 
and  the  classic  beamfonning  spatial  spectrum  is  given  by: 
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p  _aH(9,0)Ra(9,0) 

BF  a  H  (9, 4>)a(9,  (f>) 

where  R  is  the  estimation  of  the  sample  covariance  matrix. 


8.2.2  Constrained  Beamformer  (Capon  Beamformer). 

The  conventional  beamfonner  has  the  limitation  that  it  cannot  resolve  two  sources 
spaced  closer  than  the  beamwidth  of  the  array.  Modifications  to  improve  this  limitation  have 
been  proposed  by  many  researchers.  One  of  the  best  known  is  the  constrained  (or  minimum 
variance  distortionless  response  -MVDR)  beamfonner  proposed  by  Capon  [4].  The 
optimization  problem  for  the  output  power  is: 

argminP(w)  =  argminE|wHx(t)  r 

W  W  r  J 

such  that :  wlIa(9,(|))  =  1 

The  optimal  w  can  be  obtained  by  using  the  Lagrange  multiplier  technique: 

MVDR  h  /a  -1  /A 

a  (9,4>)R  a(9,<W 

and  the  Capon  (MVDR)  beamfonner  beam  pattern  is  given  by: 

Pmvdr  §)  =  h/a  j-yd-i  To 

a  (9,  (f))R  a(9,  cf)) 


8.3  Subspace  DO  A  Estimation  Methods 

Subspace,  or  eigen-structure,  methods  were  developed  from  the  spectral  decomposition 
of  the  covariance  matrix,  the  so-called  Pisarenko’s  hannonic  decomposition  method,  and  were 
first  introduced  via  the  MUSIC  algorithm  [5].  Many  subspace,  or  eigen-structure,  methods  are 
proposed  in  the  literature  due  to  its  high-resolution  property  [6]  [7]  and  the  MUSIC  algorithm 
has  been  received  more  attention. 

8.3.1  Subspace  Properties. 

If  d  signals  impinge  on  an  antenna  anay  comprising  M  sensors,  then  the  array  output 
covariance  matrix, 

Rxx  =E{XXH}=  A(9,b)PAH(9,b)  +  o2I 

where  A(9,4>)  is  the  M  x  d  steering  matrix,  P  is  the  d  x  d  signal  covariance  matrix. 

Let  d’  denote  the  rank  of  signal  covariance  matrix  P.  Then,  if  d’<  d  some  of  the  signals 
are  coherent,  (i.e.,  one  signal  is  a  copy  of  another  with  a  complex  constant,  such  as  multipath 
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signals)  and  the  signal  covariance  matrix  P  will  not  be  full  rank,  d,  but  rather  with  rank  d’.  Let 
the  eigen-decomposition  of  Rxx  be: 

M 

k=l 

where  the  eigenvalues  are  ordered  as  A-i  >  X2  >  . .  .>  AM  and  are  orthogonal  and  e*kei  =  8 
(k,l)  due  to  the  Rxx  being  Hennitean.  The  rank  of  A(0,(|))P  A(0,c())  is  d’,  and  its  null  space  has 
M-d’  eigenvalues  equal  to  &  with  the  remainder  d’  eigenvalues  larger  than  noise  variance. 
Then,  Rxx  can  be  partitioned  as: 

RXX=ESASES*+  a2  EnEn* 

where  As  =  diag(A, ,A2,  . .  .,Xd  )  and  the  corresponding  space  spanned  by  Es  is  called  the 
signal  subspace  and  En  is  the  noise  subspace. 

Subspace  methods  are  based  on  the  following  observations: 

R(ES)  c  R(A(0,4>))  and  Es  1  En 

where  R(  )  denotes  the  range  of  space  and  A  denotes  “  orthogonal  to”. 


8.3.2  MUSIC  (Multiple  Signal  Classification) 

MUSIC  (Multiple  Signal  Classification)  was  one  of  the  first  subspace  methods  to  be 
introduced  to  the  DOA  estimation  problem.  In  the  MUSIC  method,  it  is  assumed  that  the 
signal  covariance  matrix,  P,  is  full  rank,  which  means  d’  =  d  and  that  there  are  no  coherent 
signals  coming  to  the  array.  Due  to  R(ES)  cz  R(A(0,()))  and  Es  A  En ,  we  have 

a*(0k,c()k)  E„  =  0  k=l,  ...  ,  d 

and  the  noise  space  En  is  obtained  by  choosing  the  eigenvectors  corresponding  to  the 
(m-d)  smallest  eigenvalues  of  Rxx.  The  MUSIC  “spatial  spectrum”  is  then  defined  as: 

P  (6  4,)=  a"(9.tia(e,») 

uus,a  ■*’  a"(0,#.EWe,») 

The  MUSIC  estimates  of  the  DO  As  are  obtained  from  the  d  largest  peaks  in  the  spatial 
spectrum  with  searching  in  space  |0,(|)}.  When  d’<  d,  (i.e.,  there  are  coherent  signals)  the  true 
steering  vectors  are  not  orthogonal  to  the  noise  subspace.  The  MUSIC  algorithm  cannot 
handle  this  scenario  and  will  not  provide  correct  DOAs. 


8.3.3  ESPRIT  (Estimation  of  Signals  Parameters  via  Rotation  Invariance  Techniques) 

ESPRIT  (Estimation  of  Signal  Parameters  via  Rotational  Invariance  Techniques)  [7] 
is  a  computationally  efficient  and  robust  method  of  DOA  estimation.  This  method  is  based  on 
a  particular  type  of  array  geometry,  which  uses  two  identical  arrays,  see  Fig.  8.2,  to  fonn 
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matched  pairs  with  an  identical  displacement  vector,  i.e.,  the  second  element  of  each  pair  ought 
to  be  displaced  by  the  same  distance  and  direction  relative  to  the  first  element.  If  each 
subarray  has  M  sensors  it  is  not  necessary  for  the  whole  array  to  have  2M  sensors  as  sensors 
may  be  common  between  the  two  subarrays.  If  two  elements  are  separated  by  a  displacement 
d0  and  the  signals  induced  on  the  ith  pair  are  denoted  by  x,(t)  and  y,(t),  then  we  have: 

yi(t)  =  xi(t)ej^ 

where  |li  =  27tdo£Ac  and  g  =  sin0,c"JO1 ,  a  measure  of  the  DOA  of  the  signal,  and: 

x(t)  =  A(e,4>>(t)  +  nx(t) 
y(t)  =  A(e,4>)®s(t)+ny(t) 

where  O  is  an  MxM  diagonal  matrix,  with  its  mth  diagonal  element  given  by: 

<I>  =  e-^"1 

mm 

and  A(0,(|))  is  the  M  x  N  steering  matrix  (with  N  steering  vectors  corresponding  to  the  N 
directional  sources),  s(t)  denotes  the  N  source  signals  induced  on  a  reference  antenna  element, 
and  nx(t)  and  ny(t)  denote  the  noise  induced  on  the  elements  of  two  subarrays. 


Fig.  8.2  Illustration  of  ESPRIT  array  geometry 

With  eigen-decomposition  applied  to  the  correlation  matrices  Rxx  and  Ryy,  two  MxN 
matrices,  Exs  and  Eys  with  their  columns  denoting  the  N  eigenvectors  corresponding  to  the 
largest  eigenvalues  in  the  matrices  Rxx  and  Ryy,  are  obtained.  They  span  the  same  N- 
dimensional  signal  subspace  and  are  related  by: 

E  =  'PE 

ys  xs 

where  'P  is  a  unique  nonsingular  transformation  matrix.  Similarly,  matrices  A(0,4>)  and 
A(0,4>)<I>  are  related  by  another  unique  nonsingular  matrix  T,  as  the  same  signal  subspace  is 
spanned  by  these  steering  vectors.  So, 

Exs=  A(0,<|>)T  and  Eys  =  A(0,^))OT 

If  the  signal  steering  matrix  has  full  rank,  i.e.,  rank{A(0,(|))}=N,  we  have: 
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m  1  =  o 

which  implies  that  the  eigenvalues  of  ¥  are  equal  to  the  diagonal  element  of  the  <D 
matrix  and  that  the  columns  of  T  are  eigenvectors  of  ¥.  An  eigen-decompostion  of  ¥  gives 
its  eigenvalues,  and  equating  them  to  ®  leads  to  a  DOA  estimation.  If  An  is  an  eigenvalue  of 
¥,  then: 


^  =  sin 


Arg(^n) 
2  red .. 


n  =  1,2, '"jN 


8.3.4  WSF  (Weighted  Subspace  Fitting) 

The  Weighted  Subspace  Fitting  (WSF)  approach  [6]  is  a  Maximum  Likelihood  (ML)  - 
like  algorithm  based  on  subspace  decomposition  which  does  not  use  the  orthogonality  between 

2 

min  Es-  A(0,cb)T 

e,<t>,T 

W 

the  noise  subspace  and  the  steering  vector  directly.  Instead,  it  tries  to  fit  an  estimate  of  the 
signal  subspace  to  the  parameters  that  are  of  interest  using  a  ML-like  minimization.  Due  to 
R(ES)  c:  R(A(0)),  the  WSF  approach  relies  on  the  following  criterion,  which  gives  the  best 
weighed  least  squares  fit  of  the  subspaces  Es  and  A(0,4>): 

where  ||A||  B,  denotes  Trace  {AWA*},  and  T  is  an  arbitrary  (dxd’)  matrix.  The 
weighting  W  is  a  (dxd’)  Hennitean  positive  definite  matrix.  It  is  possible  to  explicitly  solve 
the  above  equation  with  respect  to  T  and  the  solution  is  given  by: 

T  =  (A*  A)-1  A*Es  =  A+Es 

and 

A  2  A  2 

V  (0,  cf))  =  ES-AA+ES  =  n^Es  =  Tr{IIAEsWEs*} 

w  w 

where  II  j  =  I-  IIA  =1-  AA+  is  the  orthogonal  projector  onto  the  null-space  of  AH. 

The  estimate  of  (0,cf))  is  obtained  as  the  minimizing  argument  of  V(0,<j)),  i.e: 

(0,  <j>)  =  arg  min  V  (0,  (f>) 

0,(J) 

Different  choices  of  the  weighing  matrix  W  lead  to  a  whole  class  of  estimates.  The 
optimal  weighting  matrix  Wopt  can  be  given  by 


Wopt=(As-a2I)2A;' 
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A 

The  unknown  quantities,  As  and  <x,  can  be  replaced  by  their  estimation,  As  and  G2 
respectively,  without  effecting  the  asymptotic  properties  of  estimate. 


8.4  Maximum  Likelihood  Estimation  Method  (MLEM) 

The  Maximum  Likelihood  Estimation  Method  (MLEM)  was  introduced  by  Fisher  [8] 
as  a  general  estimator  in  statistics  and  has  now  become  a  powerful  estimator  in  the  signal 
processing  area.  ML  estimates  have  strong  statistical  properties  and  the  ML  estimator  can 
attain  the  Cramer-Rao  Lower  Bound  (CRLB)  asymptotically  with  the  number  of  samples.  The 
CRLB  is  a  lower  bound  on  the  estimation  error  variance  for  any  unbiased  estimator. 

From  the  array  data  model  we  made  before,  the  array  response  data,  X(t),  given  J 
samples,  becomes  a  stationary,  zero-mean,  complex  Gaussian  process  with  unknown 
parameters  (0,0),  the  direction  of  arrival  of  the  signal,  and  P  the  signal  covariance  matrix.  The 
probability  density  function  is  given  by: 

j  -M  i 

p(X(i),  •  •  • ,  X(  J)  /(0, 4>),  S(t),  O2 )  =  n  Ka2  exp( — r|x(t)  -  AS(t)f ) 

t=i  ° 

where  (0,(f>)  is  the  directional  information,  S(t)  is  the  transmitted  signal  and  o2  is  the 
variance  of  the  noise  process.  The  ML  estimates  of  these  unknowns  are  calculated  as  the 
maximising  arguments  of  p(X(t)/  (0,0), S(t),  o2),  the  rationale  being  that  these  values  make  the 
probability  of  the  observations  as  large  as  possible.  Alternatively,  it  is  possible  to  minimize 
the  negative  log-likelihood  function  that  is  given  by: 

-  ln(p(X(t)  /(0, 0),  S(t),G2 ))  =  M  In  a2  +  4r  E  ||x(t)  "  AS(t)f 

G  J  t=1 

Obviously,  the  estimate  for  the  signal  waveform  is: 
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The  following  non-linear  optimization  problem  is  then  obtained  as  an  estimator  for 

(e,4>): 

(0,<j))  =  argminTrlP^A  Rxx} 

0,4> 

Maximum  Likelihood  Estimation  is  a  parametric  method  and  hence  its  resolution  is  not 
limited,  as  is  the  case  for  the  conventional  beamformer.  However,  a  multidimensional  search  is 
required  to  find  the  estimates,  resulting  in  a  high  computational  complexity.  The  ML  estimator 
presented  here  can  be  classified  as  a  deterministic  ML  estimator,  because  the  impinging 
multipath  rays  of  both  the  desired  signal  and  the  interferers  are  modelled  deterministically.  It  is 
also  possible  to  model  the  interfering  sources  as  coloured  Gaussian  noise.  For  the  detenninistic 
model,  the  number  of  signal  waveform  parameters  grows  as  the  number  of  samples  increases, 
implying  that  they  cannot  be  estimated  consistently. 

8.5  Conclusion. 

In  this  section  we  have  reviewed  some  of  the  main  estimation  methods  that  may  be  used 
to  simultaneously  detennine  the  DOAs  of  a  number  of  signals.  In  practice,  however, 
beamforming  is  of  limited  interest  in  DOA  estimation  as  its  resolution  is  limited  to  that  of  the 
antenna  array  structure.  We  will,  therefore,  for  the  remainder  of  this  study  concentrate  on  the 
use  of  Subspace  and  Maximum  Likelihood  estimation  methods  to  obtain  the  required  DOAs. 
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9.  Sensitivity  Analyses  of  Parameter  Estimation  Methods 

9. 1  Effects  of  Model  Errors  and  Sensitivity  Analysis  on  Subspace  Methods [26] 

For  the  general  model  errors,  the  array  response  model  can  be  represented  as: 

X(j)  =  C  ■  r  ■  A(0)  ■  S(j)  +  N(j) 

where  C  is  the  mutual  coupling  matrix,  which  models  the  effect  of  mutual  coupling 
between  antenna  sensors  within  the  array.  T  denotes  a  complex  diagonal  matrix  which 
includes  the  channel  gain  and  phase  errors  such  that  r„  =  a„  +  j  (ki¬ 
ln  subspace  DOA  estimation  methods,  the  array  output  covariance  matrix  Rxx  can  be 
partitioned  into  two  subspaces,  the  signal  subspace  Es  and  the  noise  subspace  En,  using 
eigenvalue  decomposition,  such  that: 

RXX=ESASES*+  a2  EnEn* 

and  R(ES)  c  R(A(0)),  which  means  that  the  signal  subspace  Es  is  spanned  by  the  array 
manifold  A(0). 

When  we  consider  array  response  errors,  the  array  manifold  becomes  A’(0),  and 

A’(0)=C  ■  T  ■  A(0)  =  C  ■  a  ■  9  ■  A(0) 

where  A(0)  is  the  true  array  manifold;  C  is  the  mutual  coupling  matrix;  a  and  (])  are 
diagonal  matrices,  the  elements  of  which  represent  the  gain  and  phase  of  each  channel 
respectively.  Therefore,  the  array  manifold  A(0)  will  affect  the  orthogonality  between  the  two 
subspaces,  Es  and  En.  In  particular,  due  to  the  effects  of  the  array  model  response  errors,  the 
real  array  manifold  A’(0)  is  different  from  the  ideal  array  manifold  A(0).  Thus  the  accuracy 
and  resolution  of  the  DOA  estimation  will  be  degraded  by  the  errors  in  the  array  manifold.  We 
refer  to  the  difference  between  the  ideal  and  real  array  parameters  as  a  model  error.  For  the 
MUSIC  algorithm,  if  vector  y  denotes  the  real  model  parameters,  e.g.,  channel  gains  a=  [ai,  0C2, 

,  aM]  or  phase  (|)  =  [()],  A . (|)m]  and  Yo  represents  the  nominal  (ideal)  parameters,  the  spatial 

spectrum  of  MUSIC  and  DOA  estimation  will  be: 

E(0;  y)  =  — - - - - - - 

AH(0;y)En(y)En(Y)HA(0;y) 

without  model  errors,  i.e.,  y  =  y0 ,  we  can  get  the  exact  estimate  of  the  true  DOAs. 

When  y^yo,  the  peaks  of  E(0;y)  will  no  longer  be  the  true  DOAs  and  the  accuracy  of  DOA 
estimations  can  be  degraded  by  even  small  model  errors.  We  will,  therefore,  consider  the 
effect  of  model  errors  on  both  accuracy  and  resolution  of  DOA  estimations. 
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9.1.1  Effects  of  Model  Errors  on  the  Accuracy  of  DOA  Estimation. 

We  assume  that  the  array  correlation  covariance  matrix  Rxx  is  exactly  estimated,  which 
means  having  an  infinite  (or  very  large)  number  of  independent  samples  (snapshots),  and  the 
effect  of  model  errors  is  sufficiently  small  so  that  the  MUSIC  spectrum  has  distinct  peaks. 
Considering  only  the  effect  of  variations  in  real  parameters,  we  can  see  how  much  the  peaks 
shift  due  to  model  error  effects.  Let: 


F(9;y)  ) — -  =  AH(6;y0)E  (y0)E  (y0)H  A(9;y0) 

E(0;yo) 

and 

F|(0;T)=« 

1  ae 

In  the  presence  of  two  sources  with  true  DOAs,  0i  and  02,  we  perfonn  the  first-order 
Taylor  expansion  to  see  how  much  the  peak  locations  shift  due  to  model  errors  and  get  the 
angular  deviation: 


3F1(0i;7)i 


A  dy  |T=To 

A0I=0_ei=  dF^fV, 


'=Yn  (y_  To) 


i  =  l,  2 

where  y-yo  is  the  model  error.  Let  us  assume  that: 


Y-Yo=°Yfi 

where  |l  is  a  random  vector  with  zero  mean  and  unit  variance  and  aY  is  a  positive 
scalar.  Thus,  we  can  show  that: 

E{(Y-Yi)}  =  0 

and  the  standard  deviation: 


std{(0  0j )}  = 


where  INI  is  the  vector  nonn. 


^(0,;^ 

dy 

y=y0 

5F1(01;y) 

30 

6=6; 

Defined,. and  <JDOA  as  follows: 
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_  Std{(0  —  0;)}  _ 

5f,(0,;y) 

5y 

Y=Y„ 

5F,(0,;y) 

50 

e  =  0; 

M  n 

6s* 


where  oDoa  represents  the  mean  square  value  of  DOA  errors  caused  by  modeling  errors 
and  N  is  the  number  of  signals.  We  can  see  that  the  angular  standard  deviation  will  be 
increased  with  Oi. 


9.1.2  Effects  of  Model  Errors  on  Angular  Resolution. 

When  the  model  errors  are  sufficiently  large,  the  algorithm  will  fail  to  resolve  two  or 
more  close  sources,  i.e.  the  model  errors  will  affect  the  DOA  resolution.  In  this  case,  as  the 
model  errors  increase  the  valley  between  two  peaks  will  become  flattened  and  finally  the  two 
peaks  merge  into  one  so  that  the  algorithm  loses  angular  resolution. 

The  valley  bottom  between  the  two  peaks  given  by  E(0,y)  can  be  seen  as  the 
minimum  of  E(0,y)  and  thus  we  have: 

F1(0*;y)  =  O, 

50 

In  general,  the  location  of  the  valley  bottom,  0  ,  will  change  with  model  error,  y.  This 
fact  complicates  the  analysis  considerably.  To  simplify  the  analysis,  we  assume  that  0*  is 

s|s 

constant  and  that  0  =0O,  where  0o=(0i+02)/2.  This  enables  us  to  study  F2(0o,y)  rather  than 
F2(0*(y),y).  Unfortunately,  this  is  only  true  for  uniform  circular  arrays  and  not  for  unifonn 
linear  arrays  in  the  presence  of  phase  errors.  Thus  the  result  of  the  first  order  analysis  here 
should  be  used  for  unifonn  circular  arrays  only. 

When  0  =0o  and  using  the  first-order  Taylor  expansion  around  yo,  the  smallest  norm 
model  error  which  will  cause  the  MUSIC  algorithm  to  fail  is  given  by: 


Y-Yo 


±F;(91i.T»)3F;^i,T>|y  =  T. 


|3F2(0„,y) 


Y  =  Yo 


and  the  failure  threshold  of  the  MUSIC  algorithm  is: 


~Vo)\\ 


where  Mis  the  number  of  the  sensors. 
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9.1.3  Conclusions  from  Sensitivity  Analyses. 

The  following  conclusions  may  be  drawn  from  the  sensitivity  analyses  with  regard  to 
the  sensitivity  parameter,  Gdoa,  and  the  failure  threshold,  Gfan: 

1)  Odoa  as  a  function  of  the  source  separation  - 

•  As  the  number  of  array  elements  is  increased,  the  sensitivity  parameter  Odoa  will  be 
decreased. 

•  The  sensitivity  parameters  for  circular  array  gain,  ogainDOA,  and  phase,  ophaseD0A,  will 
decrease  almost  linearly  as  the  source  separation  A  increases,  i.e.  ogainooA  ~  l/A,GphaseDOA~l/A. 

•  The  sensitivity  parameter  for  linear  array  gain,  ogamooA,  will  decrease  almost  in  the 
same  way  as  in  the  circular  array.  But  the  sensitivity  parameter  for  phase  is  very  small 
(ophaseD0A«  1)  and  will  keep  constant  as  the  source  separation  A  increases. 

•  For  two  adjacent  signals,  the  sensitivity  parameters  for  a  linear  array  are  smaller  than 
those  for  a  circular  array. 


2)  Gdoa  as  a  function  of  the  element  spacing  (array  aperture)  - 

•The  sensitivity  parameter  for  linear  array  phase, ophaseDOA,  is  proportional  to  1/d,  where 
d  is  the  element  spacing,  and  also  ophaseD0A  ~1/A.  So  we  have  GphaseDoA  ~l/(dA),  where  d  and 
A  are  the  element  spacing  and  source  separation  respectively.  Let  8  be  the  ratio  of  the  source 
separation  in  degrees  and  the  beamwidth,  i.e.  8  ~  (dA),  then: 

GphaSeD0A~l/5 

this  equation  shows  that  if  the  separation  between  two  sources  is  same  in  terms  of 
beamwidth,  a  linear  array  will  have  the  same  phase  sensitivity  as  a  circular  array  and  its  gain 
sensitivity,  GgainDoA,  will  decreased  as  d  increases  in  proportion  tol/d-  k  (  2<  k  <2. 5). 

•  The  sensitivity  parameters  of  a  circular  array,  GgamDOA  and  GphaseDOA  will  decrease 
rapidly  as  the  element  spacing  d  increases,  in  proportion  to  1/d-  k  (2<  k  <2. 5)  also. 


3)  Gfaii  as  a  function  of  source  separation  - 

•  The  ways  the  peaks  shift  because  of  phase  errors  are  different  in  linear  and  non-linear 
arrays.  In  a  linear  array  phase  errors  cause  the  peaks  to  shift  in  the  same  direction  by  almost 
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the  same  amount  whereas,  in  a  non-linear  array  the  peaks  move  towards  each  other  until  they 
eventually  merge  into  one. 

•  The  failure  threshold  Ofaii  is  approximately  proportional  to  A  ,  i.e.,  Ofaii  ~  A  . 

4)  Ofaii  as  a  function  of  element  spacing  - 

•  The  gain  failure  threshold  ogainfaii  is  approximately  proportional  to  d2, 

i.e.,  ogainfaii  ~  d2  and  for  a  circular  array,  the  phase  failure  threshold  is  given  by 
<7phasefaii  ~  dk  (2<  k  <  3) 

•  cjgamfaii  ~  82  and  ophasefaii  ~  82  dk~2  (circular  array),  where  8  is  ratio  of  the  source 
separation  to  beamwidth. 

Summarizing  the  results  above,  we  get  the  following  general  conclusions: 

Linear  Array:  It  has  a  low  sensitivity  to  phase  errors  and  a  larger  sensitivity  to  gain 
errors,  similar  to  a  circular  array.  When  the  source  separation  is  small  the  sensitivity  to  phase 
error  is  constant. 

Circular  Array:  The  sensitivity  to  both  gain  and  phase  errors  is  inversely  proportional 
to  the  source  separation  and  the  failure  threshold  is  approximately  proportional  to  the  square  of 
the  source  separation. 

As  the  array  aperture  (the  number  of  elements)  increases,  the  sensitivity  of  DOA 
accuracy  to  model  errors  will  decrease  and  the  failure  threshold  will  increase. 

These  conclusions  are  only  valid  for  the  assumption  that  the  data  is  infinite.  In  the 
finite  data  case,  MUSIC  is  expected  to  fail  in  the  presence  of  smaller  model  errors.  If  the  SNR 
or  the  number  of  snapshots  is  too  small,  MUSIC  will  fail  even  in  the  absence  of  model  errors. 


9.2  Effects  of  Model  Errors  and  Sensitivity  Analysis  on  Maximum  Likelihood 
Methods 

The  maximum  likelihood  DOA  estimation  is: 

0  =  arg  max  Tr  {PA(0)  R  xx } 

0 

where  Pa(9)=A(0)(Ah(0)A(0))'1  Ah(0)  is  a  projection  operator  and  is  spanned  by  the 
columns  of  A(0).  The  maximum  likelihood  estimation  of  0  involves  multiple  dimensional 
searching,  i.e.  in  the  N  dimensional  parameter  space,  if  there  are  N  signals  coming  to  the  array. 
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Here  we  will  discuss  how  the  ML  estimate  perfonnance  is  affected  by  the  DOA 
deviation,  i.e.  the  error  between  estimates  of  the  DOAs  and  the  true  DOAs  caused  by  model 


errors. 


9.2.1  The  Sensitivity  of  DOA  Estimate  Accuracy  to  Model  Errors. 

In  the  presence  of  model  errors,  the  array  output  covariance  matrix  is: 

Rxx(Y)  =  A(0;y)RsAh(0;Y)  +  o2I 

where  A(0;y)=[A(0i;y),  A(0o;y),  •••,  A(0n;y)L  and  Y represents  the  error  parameter,  i.e., 
gain,  phase  or  sensor  location  errors. 

The  ML  method  is  to  search  for  the  maximum  of  the  likelihood  function: 

L(0,Y)  =  tr(PA(e;TolRxx(Y)) 

where  PA(e;7o)=A(0;Yo)(A(0;Yo)HA(0;Yo))'1  A(0;yo)H  When  y0  =  y,  the  peaks  of 
maximum  likelihood  function  will  correspond  to  the  true  DOAs.  Otherwise,  if  Yo  ^  Yo  ,  the 
peaks  will  be  shifted  from  the  true  DOAs. 

Let  the  first-order  derivative  of  the  maximum  likelihood  function  equal  zero,  we  can 
then  get  the  maximum  of  the  maximum  likelihood  function.  That  is,  let: 


M0;Y): 


3L(0;y) 


then,  setting  Li(0;Yo)  =  0,  we  get  the  following  relationship  for  the  variation  in  the 


DOA  estimate  versus  the  model  errors: 


0-0o  =-( 


3Li(0;Y)i  dL^Y), 


,(Y-Y0) 


0  v  00  I  «o .Yo  ^  0y  |  Ho>7o  v  *  'U' 

where  (0  -  0o)  is  the  deviation  of  the  estimate  of  the  DOA  and  (j-jo)  is  the  model  error. 
The  nonn  of  the  DOA  errors  is  given  by: 


0-eo "  =(Y-Yo)m(Y-Yo) 


where 


dLjffly)  T  T  3L,(0,y) 

"  ^  dy  30  0o-7°;  1  30 


-i  3Li(0,y)i 


Let  us  assume  that  y-yo  =  oY|i,  where  |i  is  a  random  vector  with  zero  mean,  unit 
variance  and  oY  is  a  positive  scalar.  Then 


cov{0-0o}  =  o:( 


23L1(0.y),  w  3^(07)  3L1(0.y), 


3MRY) 

1  30  e°'T°; 
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Using  the  fact  that  tr{AB}=  tr{BA},  it  follows  that: 

tr{cov(0-0o)}  =  o^tr{M} 


Note  that  tr{cov(0-0o)}/N  is  the  average  variance  of  the  DOA  errors,  i.e. 


~7  tr{cov(0  -  0O )}  =  ^7  S  var(0  “  9o ) 


Therefore,  we  have  the  DOA  sensitivity  of  the  ML  algorithm  to  model  errors, 

|tr{cov(0-0o)}  _  / tr {M} 


and  this  angular  deviation  will  increase  as  tr{M}  increases. 


9.2.2  Failure  Threshold  of  ML  Algorithm. 

As  the  model  error  increases,  the  ML  algorithm  will  fail  to  provide  meaningful  DOA 
estimates  if  the  peaks  of  the  maximum  likelihood  function  occur  on  the  failure  line 
characterized  by: 

0  =  al 

where  I  =[1,1]T  and  a  is  an  arbitrary  scalar.  We  can  then  obtain: 


al-0o  =-( 


9^(0, 7)  I 


-i  3Li(0,7)| 


50  1 60  To  5y 


(Y-Yo) 


This  equation  can  be  resolved  non-uniquely  for  y  -  yo,  given  any  value  of  a.  We 


perfonn  a  singular  value  decomposition  of 


3L,  (0,y) 


,  to  give: 


dLi(0,Y)i 


=  U,2IV1H 


Define  8  =  V/1  (0  —  0O)  and  then: 


al-0o  = 


3L,(e,y) 


U.2,8 


or: 
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31^9,7) 

ae 


UiZjI 


The  minimum  norm  solution  of  this  under-determined  set  of  equations  could  be 
obtained  by  various  matrix-solving  routines.  We  only  need  to  find  the  smallest  ||8||  and 
smallest  ||y-yo||-  Due  to  V i  being  an  orthogonal  matrix,  the  nonn  of  8  is  same  as  the  norm  of  y- 
Yo.  Thus,  if  we  consider  the  effect  of  model  errors  on  each  element,  the  failure  threshold  of 
the  ML  algorithm  will  be  related  to  ||8||  and  we  have  the  failure  threshold  of  the  ML  algorithm 
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10.  Selection  of  Estimation  Method  and  Algorithm  Development. 

10.1  Selection  of  Estimation  Method. 

In  section  8  of  this  report  we  reviewed  the  main  methods  of  parameter  estimation; 
beamforming,  subspace-based  methods  and  maximum-likelihood  methods.  In  practice, 
however,  beamforming  is  of  limited  interest  in  DOA  estimation  as  its  resolution  is  limited  to 
that  of  the  antenna  array  structure.  In  section  9  of  this  report  sensitivity  analyses  were 
produced  to  investigate  the  effect  of  model  errors  on  both  Subspace  methods  and  Maximum 
Likelihood  methods.  As  a  result  of  these  analyses,  it  was  decided  that  in  the  case  of  known 
frequency  an  ML  method  would  be  considered  and  for  unknown  frequencies  subspace  methods 
would  be  investigated,  as  these  would  give  the  best  angular  resolution  while  having  a  lower 
computational  complexity  than  ML  methods.  A  number  of  estimation  algorithms  using  sub¬ 
space  methods  were  considered,  with  special  emphasis  on  the  solution  of  multi-parameter 
estimations,  and  details  follow  in  this  section. 

10.2  Algorithm  Development. 

In  section  8.3  of  this  report  we  reviewed  the  basics  of  Subspace-based  estimation 
methods  with  reference  to  the  MUSIC,  ESPRIT  and  WSF  methods  in  particular.  As  part  of 
this  project  work  has  been  carried  out  into  using  various  subspace  methods  to  estimate  a 
number  of  parameters  simultaneously.  Maximum  Likelihood  methods  were  considered  in 
Section  8.4. 

10.2.1  Maximum  Likelihood  Method. 

We  will  consider  an  ML  method  to  estimate  the  Mutual  Coupling  Matrix  in  the  case  of 
known  signal  frequency.  For  the  papers  of  mutual  coupling  modeling  and  previous  work,  see 
[28],[30],[38],[43],[50],[51],[53],and[58], 

We  make  the  standard  assumptions  underlying  the  algorithms: 

(1)  The  signal  and  noise  processes  are  stationary  and  ergodic  over  the  observation 

period. 

(2)  The  columns  of  A(9)  are  linearly  independent. 

(3)  The  noise  is  uncorrelated  with  the  signal,  and  it  has  zero  mean  value.  Its 
covariance  matrix  is  full  rank  and  is  known  except  for  a  multiplication  constant  o“. 

(4)  The  total  number  of  sources  is  known  (2  <  p  <  M  )  and  the  sources  have  the  same 
known  frequency. 
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(5)  The  MCM  satisfies  rank(C)  =  M  which  implies  that  rank  (CA(0))  =  rank(A(0))  =  p 

10.2.1.1  Estimation  process 

The  algorithm  proposed  here  consists  of  two  estimation  procedures:  estimation  of  signal 
directions  of  arrival  {0n}  and  estimation  of  mutual  coupling  matrix  C. 

1.  Estimation  of  directions  of  arrival  {0n} 

Defining: 

Q,  =E||x(j)-c°A(e°)s«(j) 2 

j=l 

S„(j)=(BjB„)''B0“x(j)  ...(10.1) 

where  ||  •  ||  denotes  the  Euclidean  nonn.  J  is  the  number  of  samples  (snapshots).  Bo  = 
C°A(0u)  denotes  the  array  steering  matrix  with  the  effects  of  mutual  coupling,  C°  and  0°  are  the 
initial  values  of  the  mutual  coupling  matrix  and  DOAs  respectively. 

During  this  first  step,  the  algorithm  perfonns  successive  minimization  operations  on 

A 

each  column  of  B0  and  holds  all  of  other  columns  and  associated  components  of  So(j)  fixed. 
Suppose  that  we  perfonn  the  minimization  with  respect  to  the  kth  column  vector,  Qi  can  be 
written  as: 

Q,  =L||xk(j)-bo(el)s„t(j):  ....(io.2) 

j=i 

h  A  h  A  k 

where  bo(0k)  is  the  kth  column  of  B0 ,  Sor  (j)  is  the  kth  component  of  So  (j) ,  and  Xk(j)  is 
given  by: 

Xk  (j)  =  X(j)- c°A(e»  )sS  (j)  =  X(j)- B0S5  (j)  ....(10.3) 

A 

where  Sk(j)  is  a  simplification  of  So(j)  with  the  k  part  replaced  by  zero. 

The  minimization  of  (10.3)  with  respect  to  ao(0k),  using  (10.1)  with  B0  replaced  by 
bo(0k),  is: 

ek  =argmin£  X k (j)- b0 (ek )(b0 (o t )H b0 (e„ )f 1  b0 (e„ )H Xk (j)| ' 

0k  j=l 

which  is  also  equivalent  to: 

0k  =argmax-i-^|xk(j)Hbo(0k)  ....(10.4) 

yk  JYl  :=1  1 
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where  M  =  bo(0k)Hbo(0k) 

To  maximize  (10.4),  we  can  search  over  the  space  of  {0k}  and  finally  get  the  directions 
of  arrival  (DO  As)  {0n}.  Based  on  these  DO  As,  we  can  reconstruct  the  new  steering  matrix 
A(O')  and  put: 

S»,0)=(Bo“Bo,t'Bo"1X(j) 

where  Boi  =  C°  A(0' ),  01  =  [0i  02  ■  •  •  0P  ] .  The  initial  mutual  coupling  matrix  C°  of  this 
first  step  is  fixed  and  we  have  to  iterate  and  update  C  at  the  second  step  to  approach  the  true 
DOAs  and  MCM  C. 


2.  Estimation  of  the  mutual  coupling  matrix  C. 

At  this  step,  we  can  estimate  the  MCM  with  the  DOAs{0n}provided  by  the  last  step. 
Since  we  have  new  values  for  Soi(j)  and  A(0 1 ),  Q  can  be  written  as: 

Q  =  Z||x(j)-CA(el)So,(jf  ....(10.5) 

j=l 

Because  C  is  an  M  XM  banded  complex  symmetric  Toeplitz  matrix  in  the  case  of  a 
unifonn  linear  array  or  an  M  XM  complex  circulant  matrix  in  the  case  of  a  uniform  circular 
array  and  A(01)SoiO)  is  an  M  X 1  complex  vector,  the  following  two  lemma  [11?]  will  be 
useful. 


Lemma  10.1  For  any  M  X 1  complex  vector  x  and  any  M  XM  banded  complex 
symmetric  Toeplitz  matrix  B: 

B.x  =  Qi(x).b 

where  the  L  X 1  vector  b  is  given  by: 

bi  =  Bu,  i  =  1,2,  ...,  L 

where  L  is  the  highest  super-diagonal  that  is  different  from  zero.  The  M  XL  matrix 


Qi(x)  is  given  by  the  sum  of  the  following  matrices: 

rTl  _  Jxp+cH  P  +  q^M-1 
1  pq  [o  otherwise 


xP-q+i  p^q^2 

0  otherwise 


Lemma  10.2  For  any  M  X 1  complex  vector  x  and  any  M  XM  complex  symmetric 


circulant  matrix  A: 
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A.x  =  Q2(x).a 

where  the  component  of  the  L  X 1  vector  a  is  given  by: 

ai  =  An,  i  =  1,2,  L 

where  L  =  M/2  +lwhen  M  is  even  and  L  =  (M  +l)/2  when  M  is  odd. 

The  M  XL  matrix  Q2(x)  is  given  by  the  sum  of  the  following  four  M  XL  matrices: 

rTt  fVc-i  p  +  q  <  M  - 1 
'  pq  [o  otherwise 

[T]  j>W  P^2 

pq  [0  otherwise 

[Tj  =  -  P<C>S1 

pq  (0  otherwise 

[T  1  =  rxp+q-M-i  2<q<l,p  +  q>M  +  2 

4  pq  [0  otherwise 


[T.]pq  = 


Lp+q-M-l 


where  1  =  M/2  for  even  M  and  1  =  (M  +l)/2  for  odd  M  . 


Using  Lemma  10.1  or  Lemma  10.2,  depending  on  the  array  configuration,  we  have: 


CA(e%(j)=cw 


CW  =  T  G 

AMxL'j  Lxl 


where  the  L  X 1  vector  G  is  given  by: 

G,=Ch, 


i  =  1,  2,  •••,  L 


by  letting  CW  =  TG,  we  can  place  the  unknowns  in  the  L  X 1  vector  G  and  the  known 
values  in  the  M  XL  matrix  T.  Thus,  (10.5)  can  be  rewritten  as: 

Q  =  L||x(j)-TG||2 


and  vector  G  as: 


g=I^(tht)-|thxO) 


and  we  can  now  reconstruct  the  MCM  C  from  vector  G: 


Toeplitz(G)  ForULA 
Circulant(G)  For  UCA 
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10.2.1.2  Convergence  check 

During  both  the  estimation  steps,  we  check  the  convergence  of  (10.2)  and  (10.5). 
When  both  Qi(k_1)  -  Qi(k)  <  8  and  Q(k_1)  -  Q(k)  <  e  (e  is  a  given  threshold ),  the  algorithm 
has  converged  and  we  get  the  estimation  of  the  DOAs  and  the  MCM  C  .  Otherwise,  the 
algorithm  will  be  repeated  between  the  two  steps  until  Qi  and  Q  converge.  Note  that  the  cost 
functions,  Qi  (Qi  >  0)  and  Q  (Q  >  0),  are  convergent  series  and  the  convergence  of  both  is 
guaranteed. 


10.2.2  Iterative  MUSIC  Method 

In  [9],  an  iterative  subspace  method  based  on  the  MUSIC  algorithm  was  presented  to 
deal  with  the  estimation  of  DOAs  and  unknown  mutual  coupling  as  well  as  gain  and  phase 
errors.  The  algorithm  first  searches  for  peaks  in  the  MUSIC  spectrum,  similarly  to  the  original 
MUSIC  algorithm  without  mutual  coupling,  to  get  an  estimation  of  the  DOAs  and  then 
estimates  the  mutual  coupling.  This  procedure  is  repeated  several  times  and  becomes  a  two- 
step  iterative  optimization  for  the  cost  function  [9]: 

J.  =  Z||Ccra(eJ! 

P=i 

A 

where  E„  are  the  noise  eigenvectors,  eigen-decomposed  from  the  estimation  of  the  data 
covariance  matrix  R,  C  is  the  mutual  coupling  matrix  and  T  is  the  gain  and  phase  matrix,  in 
this  case  it  is  assumed  that  T  =  I  as  only  the  mutual  coupling  is  being  considered. 

Step  1 :  Estimating  the  DOAs. 

Search  for  the  N  highest  peaks  of  the  MUSIC  spectrum  defined  by: 

E„Ca(e|‘J 

where  the  initial  mutual  coupling  matrix  C  can  initially  be  chosen  as  an  identity  matrix 
and  the  estimations  of  the  DOAs  are  associated  with  the  N  highest  peaks  in  the  spectrum. 

Step  2:  Estimating  the  Mutual  Coupling  Matrix 

In  this  step,  the  vector  a(0P)  in  the  cost  function  Jc  will  be  updated  by  the  estimated 
DOAs,  {0p}(p  =1,  2,  . . .,  N),  and  then  the  mutual  coupling  matrix  C  will  be  estimated  by 
minimizing  the  cost  function  Jc . 


:(e)  = 
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It  can  be  shown  that  for  a  Uniform  Linear  array  (ULA)  the  MCM,  C,  is  Toeplitz  and  for 
a  Uniform  Circular  Array  (UCA)  the  MCM  is  Circulant.  The  cost  function  Jc  can  then  be 
written  as: 

Jt  =  i»(eP)HcHE„E;ca(ep) 

P=1 

=  c"|ZQ(p)"e],e,"q(p)Jc 

where: 

Q(p)  =  Q(a(0p))  andci  =  Cii,  i=l,2,...,L 

and  Q(X)  and  L  are  dependent  on  the  array  configuration. 

The  minimization  of  cost  function  Jc  is  a  quadratic  minimization  problem  under  a  linear 
constrain  where  the  MCM  C  is  assumed  to  be  unity  on  the  main  diagonal  (Cu  =1).  Thus  using 
Lagrangian  multipliers,  the  constrained  minimization  problem  can  be  solved  as: 

c  =  G'1w(wTG_1w)*1 

where  G  is  the  matrix: 

G  =  |]Q(p)HEnE;Q(p) 

P=1 

and  vector  W  =  [1  0  •••0]  represents  the  linear  constraint. 

A 

Finally,  the  MCM  C  can  be  reconstructed  from  the  estimated  vector  c  .  Note  that  the 
cost  function  Jc  will  be  decreased  at  each  iteration  and  the  convergence  of  the  algorithm  is 
guaranteed. 


10.2.3  Separable  Sub-space  Method 

A  following  method  for  the  joint  estimation  of  signal  frequencies,  DO  As  and  sensor 
mutual  coupling  was  presented  as  a  paper  at  the  IEEE  34th  Asilomar  Conference  on  Signals, 
Systems  and  Computers,  October  29  -  November  1,  2000,  in  Pacific  Grove,  California  [10]. 
This  paper  is  attached  to  this  report  as  Appendix  A. 

Consider  an  array  composed  of  M  sensors  with  each  sensor  output  being  fed  to  a  tap- 
delay-line  with  m  taps.  The  delay  between  adjacent  taps  is  t0.  Let  Xi(t),  (1  =1  ,2  ,  •  •  •  ,M  ), 
denote  the  output  of  the  1th  sensor  at  time  t  and  X(t)  =  [x(t),  •••,  x(t  -(m  -l)to)]T . 
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Assume  that  p  narrow-band  sources  with  different  frequencies  and  directions  impinge 
on  the  array  and  the  incoming  signals  can  be  divided  into  D  groups  according  to  their 
wavelength,  i.e: 

D 

p=XPk 

k=l 

where  pk  is  the  number  of  signals,  from  different  directions,  in  the  kth  group,  and  whose 
wavelength  is  A,k  (k  =1  ,2  ,  •  •  •  ,D  ). 

The  1th  sensor  (1=1  ,2  ,  •  •  »,M  )  output  vector  maybe  written  as: 

D  M 

X,  (t)  =  zz  (Cu  K  M<ok  (ek  K  (t)+ N,  (t) 

k=l  i=l 

(k=l  ,2,...,D,i=l  ,2,...,M) 

where,  bi(0k)  =  [exp(-jcOk  I,  (0k, i ))  •  •  •  exp  (-j  cok  x,  (0k,Pk))]  is  a  1  X  pk  spatial  steering 
vector;  a(o)k)  =  [  1  exp(-jo)kto)  •••  exp(-jcOk(m  -l)t0)]T is  an  m  X  1  temporal  steering 
vector;  Cy  (otk)  is  the  (l,i  )th  entry  in  the  MCM  which  represents  the  mutual  coupling  effect 
from  other  sensors  on  the  1th  sensor;  Sk(t  )  =  [Sk,i(t)  •  •  •  Sk,pk(t)]T  is  a  pk  X  1  signal  vector 
from  the  kth  group  of  narrow  band  sources  and  Ni  =  [Ni(t)  •  •  •  Ni(t  -(m  -l)to)]T  is  an  m  X  1 
additive  noise  vector  for  the  1th  sensor. 

Therefore,  the  array  data  model  can  be  described  as  follows: 

x(t)=[x,(t)T  x;(t)T  ...  XM(t)Tf 
=  £(C(c%  )B(et ))®  (a(tok  )Sk(t))+N(t) 

k=l 

where: 

'b,(ek)‘ 

B(e„)=  : 

_b„(ek)_ 

is  an  M  Xpk  spatial  steering  matrix;  N(t)  =  [  Ni(t)T  N2(t)T  •  •  •  NM(t)T]T  is  a  Mm  X 1 
additive  white  Gaussian  noise  matrix;  the  symbol  0  denotes  the  Kronecker  product  and  C(tOk) 
(k  =  1,  2,  •••,  D)  is  the  M  XM  MCM.  For  a  ULA  or  UCA,  C(ff>k)  is  either  a  Toeplitz  matrix 
or  a  circulant  matrix,  respectively  [9], 

10.2.3.1  Subspace  decompositions. 

Assume  that  the  number  of  signals  p  and  the  number  of  frequency  groups,  D,  are 


known. 
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Define: 


where: 


<=i£E[*,Mh-ooxl(t)"] 

tvt  i=i 

(h  =  1,  2,  D) 


ni  =  — jE[yi(t-(h-l)t0)Y(t-(h-l)t0)H] 
m  t~i 

(1  =  1,  2,  p) 


y,  (t  -  (h  -  l)t0 )  =  (Xj  (t))h  =  X!  (t  -  (h  -  l)'t0 ) 


Y(t  —  (h  —  l)t0 )  =  [y  j  (t  —  (h  —  l)t0 )  y2(t-(h-l)t0)  •••  yM  (t  -  (h  -  l)t0  )]T 
=  [x  i  (t  —  (h  —  l)t  o )  x2(t-(h-l)t0)  •••  xM(t-(h-l)t0)]T 

and  E  (•)  denotes  statistical  expectation.  {r’h  }°=1  and  jr|’  }f=1  are  the  frequency  and  the 
direction  vectors  respectively.  The  range  spaces  spanned  by  the  columns  of  these  vectors  are 
contained  or  are  equal  to  the  range  spaces  of  a(co)  and  B(0),  respectively.  We  have: 

)  £  ^(b(0)) 

and: 

N(r:)±R(»(<o)l  AT(n;)=«(B(0)) 

where  R  (•)  and  N  (•)  denote  the  range  space  and  null  space.  As  is  well  known, 
subspace  methods  are  based  on  the  above  geometrical  relationships.  Therefore,  we  may  use 
correlation  processing,  in  spatial  and  temporal  dimensions  respectively,  to  get  the  estimates  of 
{r’h  }°=1  and  J r)’,  and  then  compute  their  null  subspaces.  Finally,  the  frequencies,  DOAs 
and  the  MCM  can  be  estimated  with  subspace  methods  by  searching  corresponding  spaces. 


10.2.3.2  Estimation  algorithm. 


Step  1 .  Frequency  Estimation. 

(1)  Estimation  of  the  frequency  vectors  r’h  and  rh : 


K  =(l  -62eh)H 


where: 
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eh  =  0  •••  0  1  0  •••  0  , 

(h  =  1,  2,  •••,  D) 

G2  is  the  estimate  of  the  noise  variance  and  N  is  the  number  of  samples. 

(2)  Gram-Schmidt  (GS)  orthogonalization  and  fonnation  of  the  temporal  projection 
matrix  Pw: 

From  the  vector  {i\}°=1 ,  we  can  get  D  orthogonal  vectors  {qij,  (k  =  1,  . . D),  via  GS 
orthogonalization.  Let  Qto  =  [q  l  ,qe  ,  •  •  •,  qo  ],  then  compute  the  temporal  proj  ection  matrix 
Pw=  I -QtoQco  H,  which  spans  the  null  space  of  {a(tOk)},  (k=l,  ...,D). 

(3)  Estimate  the  unknown  frequencies  with  the  temporal  projection  matrix  P0) : 

The  frequencies  {<x>k} ,  (k  =  1,  . . .,  D),  are  estimated  as  the  D  largest  peaks  of  the 

function: 

p(<°)= 

a(co)  P0)a(co) 

by  searching  over  the  frequency  sector  of  interest. 

Step  2.  Direction  and  Mutual  Coupling  Estimation 

(1)  Estimation  of  direction  vectors,  rj’/  and  T|/ : 

i  m  i  N 

mtfN  tT 

f|,  =(q;-a2ei)H 

where: 

e,  =  0  •••  0  1  0  •••  0  , 

0^9  ' 

(l  =  k  2,  •••,  p) 

A 

(2)  By  GS  orthogonalization  of  {qJ^j  ,  the  p  orthogonal  vectors,  J C , }  [C ,  and  spatial 

orthogonal  projection  matrix  Pe  =  I  -  Qe  Qe  H  are  obtained,  where  Qe=[Ci  C2  ...  Cp]- 

(3)  For  each  frequency  {tOk},  (k  =  1,2,  . . .,  D),  estimated  in  Step  1,  the  directions  0; k, 
(ik=  1 ,2,  •••,  pk)  and  mutual  coupling  matrix  C(cOk)  can  be  estimated  by  using  the  Iterative 
MUSIC  method  described  in  section  10.2.2  with  the  following  equations: 
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p(e)=(B(e)“cK)“p,cK)B(e)i1 
£(“k  )  =  (gK  I  w  i''  l(',((!)  )  w)  1 

where  G(cok )  =  ^d(0;  )H  PeD(0it )  and  w  =  [1,  0,  •  •  •,  0]T .  Using  the  properties  of 

ik=1 

Toeplitz  and  Circulant  matrices  and  Lemmas  10.1  and  10.2,  matrix  D(0ik )  can  be  calculated 
[1 1]  by  letting  C(cok  )B(0lk  )=  D(0ik  >  ,  Le.,  the  unknowns  are  placed  in  vector  c  and  the  known 
values  in  the  matrix  D. 


10.2.4  Alternating  Iterative  Method. 

This  section  presents  an  alternating  iterative  based  approach  for  the  estimation  of 
DOAs,  signal  frequencies  and  corresponding  array  sensor  gain/phase  errors  from  the  received 
signals.  The  previous  work,  see  [27], [41], [44], [47], [55],  and  [59].  The  approach  consists  of 
two  parts:  first  estimate  signal  frequencies  and  then  estimate  the  directions  of  arrival  and  the 
corresponding  gain/phase  errors.  Since  this  multiple  parameter  estimation  procedure  is  a  non¬ 
linear  optimisation  problem,  the  alternating  iterativeness  will  be  used  to  get  the  solution,  i.e, 
one  can  first  estimate  frequencies  and  DOAs  based  on  supposing  that  the  array  gain/phase 
errors  are  initially  known,  which  can  be  determined  from  nominal  or  recently  measured  values, 
and  then  estimate  the  array  gain/phase  errors  with  the  estimated  frequencies  and  DOAs.  This 
procedure  will  be  repeated  until  the  cost  function  reaches  a  minimum. 


10.2.4. 1  Estimation  procedure 

1)  Estimation  of  the  signal  frequencies  (Ok ,  (k  =  1  . . .  D). 

Applying  conventional  frequency  estimation  methods,  such  as  a  periodogram  or  other 
temporal  spectral  estimation  methods  [12],  to  the  received  signal  matrix  X,  D  signal 
frequencies  can  be  estimated.  For  each  estimated  signal  frequency,  one  can  then  estimate  the 
signal  directions  of  arrival  and  the  gain/phase  matrix  with  the  following  iterative  procedures. 

2)  Estimation  of  the  DOAs  0kn  ,(n  =  1  ...  pk). 
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Suppose  that  the  initial  gain/phase  matrix  Ek'01  is  known;  this  could  be  based  on  either 
the  nominal  gain/phase  values  or  on  any  recent  calibration  infonnation.  For  the  ith  iteration, 
the  azimuth  and  elevation  space,  0k  =[0,  (|)]  ,  are  searched  with: 

DlI,(0k)  = - - - - 

|Enr^a(0k|  ...(10.6) 

k  =  1,  2,  •••,  D 

where  En  is  a  matrix  whose  columns  are  the  eigenvectors  corresponding  to  the  smallest 
eigenvalues  of  covariance  matrix  Rx  =  E  {  X  XH  } .  Then,  the  pk  peaks  indicate  the  signal 

directions  of  arrival,  {§k  )  | ,  at  the  ith  iteration. 

3)  Estimation  of  the  gain/phase  matrix  Tk. 

When  the  DOAs  with  frequency  0)k  have  been  estimated,  the  estimation  of  the 
gain/phase  matrix  at  that  frequency,  T k,  becomes  a  constraint  optimisation  problem  for  the  cost 
function: 


n=l 


k  =  1,  2,  •••,  D 

Since  the  gain/phase  error  matrix,  r\,  is  an  M  x  M  diagonal  matrix  and  a((Ok,  0kn)  is  an 
M  x  1  vector,  we  have: 

Eka(®k’®kn  )=  Q(®k’^kn  )^k 

where  8k  is  a  vector  composed  of  the  diagonal  elements  of  Tk  and  matrix  Q(ff>k,  0kn)  is 
an  M  x  M  diagonal  matrix  whose  diagonal  elements  are  the  elements  of  vector  a((Ok,  0kn). 

Therefore,  the  cost  function  Lk  can  be  rewritten  as: 

p 

K  =S»k.el.)"r1E"E.Uak,e1.) 

n=l 

=  8?|tQ(<o1,ek,)E;EnQ((o1,et-)k 

with  the  optimisation  constraint: 

5kw  =  l,  w  =  [l  0  •••  Of 

The  result  of  this  minimization  is  given  by: 
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^update)  _  Zk1w/(wTZ“1w)  ...(10.7) 

where: 

p 

A  =  £<*(“»  A>XQk.eJ 

n=l 

and  finally: 

r ( update,  =diag|g (update)  j  _(|()  8) 


10.2.4.2  Algorithm  summary. 

The  alternating  iterative  method  can  be  summarized  as: 

(1)  Estimation  of  frequencies  with  conventional  frequency  estimation  methods. 

(2)  Estimation  of  data  covariance  matrix: 


C=jZx,xI 

J  j=i 


(3)  Eigenvalue  decomposition  (EVD)  of  the  covariance  matrix  Rx  to  separated 
signal  subspace  Es  and  noise  subspace  En. 

(4)  Set  i  =  0  and  set  initial  value  of  the  gain/phase  matrix  for  each  frequency  (Ok, 


where  k  =  1,  2,  ...,  D. 


(5)  For  each  estimated  frequency  (Ok,  the  estimation  of  the  directions  of  arrival  can 
be  calculated  using  (10.6),  where  the  pk  peaks  indicate  the  directions  of  arrival,  Op,,1 11  (n=l,  . . ., 
Pk). 

(6)  Using  (10.7)  and  (10.8),  estimate  the  gain/phase  matrix  Fk(1). 

D 

(7)  Compute  the  cost  function  Lk(1)  and  the  cost  function  L  =  ^Lk.  Repeat 

k= 1 

steps  (4)  -  (7)  with  updated  DOAs  0kn<o  (  n  =  1,  . . .,pk)  and  gain/phase  matrix  Tk0  until  the 
reduction  in  cost  function  Lk(1)  -  Lk(1+1)  <  8,  (  where  8  is  a  given  threshold).  For  D  frequencies, 
this  procedure  will  be  repeated  D  times  until  the  cost  function  L  reaches  a  minimum. 


As  the  above  estimation  is  a  nonlinear  alternative  iterative  optimisation  procedure, 
whether  or  not  the  iterative  procedure  converges  to  the  global  optimum  depends  on  the  choice 
of  initial  values.  This  method  may  converge  to  a  local  optimum.  To  improve  the 
performance  of  convergence,  we  will  introduce  a  simulated  annealing  algorithm  and  then 
present  a  global  optimisation  method  for  the  estimation  of  DOAs,  frequencies  and  array 
gain/phase  error. 
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10.2.5  Global  Optimisation  Based  on  Simulated  Annealing 

The  joint  DOA  and  frequency  dependent  gain/phase  error  estimation  algorithm 
proposed  here  is  based  on  simulated  annealing  and  consists  of  two  parts.  The  first  is  to  fonn  a 
projection  matrix  from  the  covariance  matrix  and  to  then  construct  a  spectral  estimation  matrix 
and  cost  function  from  the  projection  matrix.  The  second  is  to  minimize  the  cost  function  by 
simulated  annealing. 


10.2.5. 1  The  projection  matrix  and  spectral-estimation  matrix. 

In  Section  7  we  have  described  the  modified  frequency-direction  array  response  model 
with  unknown  gain  and  phase  errors.  Our  work  here  is  first  to  fonn  an  orthogonal  projection 
matrix  P  from  the  covariance  matrix  Rx.  Such  a  projection  matrix  can  be  formed  from  the  p 
principal  eigenvectors  of  the  covariance  matrix. 

On  obtaining  the  projection  matrix  P,  we  can  construct  a  spectrum-estimation  matrix: 

S^r’  ^  ~  Tr([PA>,e)][pA>,e)H]) 

where: 

a'(o),  0) = [r(co,  )a(co,  ,  0! ):  r(co2  )a(co2  ,  e2  ):•••;  r(coD  )a(cod  ,  eD )] 

is  the  modified  frequency-direction  steering  matrix  with  the  sensor  gain/phase  errors. 

Because  the  trace  of  the  spectral-estimation  matrix  contains  infonnation  for  the  signal 
frequencies,  DOAs  and  array  gain/phase  errors,  the  estimation  of  all  these  variables  may  be 
obtained  by  minimizing  the  cost  function: 

L(r,®,0)  =  TrffPA^w,  eyfpA^o),  0)H  j) 

=  £Tr([pA'K,e)][pA'K,e)H]) 

k=l 

The  minimization  of  this  cost  function  is  a  non-linear  optimisation  problem.  We  will 
utilize  simulated  annealing  to  perform  a  global  optimisation.  Regarding  simulated  annealing 
and  its  applications,  see  papers,  e.g  [15], [18], [23], [36], [47], [49],  [57]and  [60]. 
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10.2.5.2  Global  optimisation  procedure 

As  has  been  described  above,  the  proposed  algorithm  consists  of  two  parts;  fonning  the 
cost  function  from  the  spectrum-estimation  matrix  and  then  using  simulated  annealing  to 
minimize  the  cost  function.  The  global  minimization  procedure  by  simulated  annealing  can 
be  outlined  as  follows: 

Initialization 

set  annealing  schedule  (T0 ,  n0 )  •  •  •  (Tf ,  nf )  ; 
generate  (r°,co°,90)  randomly; 
compute  L(r°,co°,00); 

where 

L(r,  to,  0)  =  £  tr{[pr(cok )  A(cok ,  0k  )][Pr(cok )  A(cok ,  0k  )]H  } 

k=l 

outer  loop 
set  i=l; 
inner  loop 

from  k=  1  to  ni  ; 

generate  random  perturb  (rk ,  cok ,  0k ) ; 

compute  L(rk ,  cok ,  0k  ) ; 

generate  random  number  c,  with  (0,1) 

uniform  distribution ; 

if  £,  <  p  or  if  A  >  0  ,  accept  (rk ,  C0k ,  0k )  ; 

where 

p  =  min(l,exp(-A/Ti)), 

A  =  L(rk -1 ,  cok-‘ ,  0k_1 )  -  L(rk ,  cok ,  0k ) , 

otherwise  ,reject; 
end  of  inner  loop; 

set  L(rk ,  cok ,  0k )  generated  by  inner  loop  to 
L(r°,co°,00); 
i  =  i  +  1  ; 

repeat  until  i  >  f  or  according  to  other  stop 
criterion  to  end  outer  loop. 
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In  the  above  procedure,  the  kernel  of  the  simulated  annealing  algorithm  is  to 
detennine  the  annealing  schedule  {(T0 ,  n 0 )  •  •  •  (Tf ,  n f )}. ;  where  T0  is  the  initial  temperature, 

T,  (i  =  1,2,  •  •  • ,  f )  is  the  current  control  temperature,  and  n  1  (i  =  1,2,  •  •  • ,  f )  is  the  cycle  number  of 
the  inner  loop. 

The  random  perturbation  (rk ,  cok ,  0k )  generated  at  the  kth  repetition  will  be  accepted 


with  probability  p: 


where: 


p  =  min(l,exp(-A/T1)) 


a  =  L(rk -1 ,  cok-‘ ,  ek_1 )  -  L(rk ,  ,  ek ) 


k  Ak ' 


L(r,  co,  0)  =  £  tr{[pr(cok )  A(cok ,  0k  )][Pr(cok )  A(cok ,  0k )]"  } 


If  A  >  0,  then  accept L(rk,cok,0k) ,  otherwise  accept  L(rk,cok,0k)  with 


k  „k  nk ' 


probability  p  =  exp(-  A/Ti ) . 

Repeating  the  above  procedure  can  reach  the  global  minimum  (or  approximately  global 
minimum)  and  the  proof  of  convergence  can  be  seen  in  [13,  14]. 

10.2.5.3  Algorithm  summary. 

The  main  procedure  of  the  proposed  algorithm  may  be  summarized  as  follows: 

1  )  Initialization  : 

Select  the  initial  values  for  the  signal  frequencies  {co0},  directions  {ou},  and  array  error 
matrix  {f0 }.  Usually  the  initial  values  can  be  selected  randomly  or  from  the  last  measured 


values. 


2  )  Use  all  the  available  data  vectors  to  estimate  the  data  covariance  matrix  : 


R.  =)jx(j)x(j)“  (j  =  1,2,  ....  ,  J 

J  j=l 


where  J  is  the  number  of  samples. 

3  )  Form  the  orthogonal  projection  matrix  P  using  eigenvalue  decomposition  (EVD) 
of  the  data  covariance  matrix  Rx . 

4)  Construct  the  spectral-estimation  matrix  S: 


S(r,co,0) 


tr([PA’(G>,0)][PA'((O,0)]H) 
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and  the  cost  function: 


L(r,co,0) 


1 

[s(r,co,e)] 


tr{[PA'  (oj,  0)][PA'  (to,  0)]H  } 


Z  tr{[pr(co,k  )A(cok ,  0)][Pr(cok ) A(cok ,  0)]H  } 


5)  Search  (r,  to,  0)  by  simulated  annealing  to  minimize  the  cost  function  and  get  the 
estimation  of  the  signal  frequencies,  directions  of  arrival  and  array  sensor  errors  (gains  and 
phases). 
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11.  Simulation  Results. 

11.1  Convergence  Properties. 

As  convergence  is  an  important  property  for  iterative  estimation  methods,  the 
convergence  rates  of  the  ML  method  of  Section  10.2.1  and  the  IMUSIC  method  of  section 
10.2.2  are  compared  in  the  first  example.  Consider  a  uniform  linear  array  with  8  sensors,  each 
sensor  separated  by  A/2,  with  2  equi-power  narrow-band  sources  of  known  frequency  located  in 
the  far  field  of  the  array  at  directions  0i  =  30°,  02  =  45°.  The  signals  are  assumed  to  be 
uncorrelated.  Additive  noise  is  injected  with  a  SNR  of  30dB  referenced  to  each  signal  source, 
i.e., 

p/a2  =  1030/10I2 

500  snapshots  of  array  data  are  accumulated.  The  initial  values  for  the  DOAs  were 
chosen  as  0|  =  35°,  02  =  50°,  and  the  initial  mutual  coupling  matrix  is  chosen  as  an  identity 
matrix.  Figure  11.1  shows  the  DOA  estimation  of  0  2  =  45°  during  35  iterations  using  the  ML 
and  IMUSIC  methods. 


Fig.  1 1.1  The  DOA  estimation  of  the  9  =45  using  the  ML  and  IMUSIC  methods  versus  the  number  of  iterations 
using  a  ULA  of  8  sensors  when  two  waves  are  incident  from  35°  and  45°  and  500  snapshots  are  taken  with  a  SNR  =30dB 

In  comparison  with  the  ML  method,  the  IMUSIC  method  has  a  slower  convergence  rate 
to  the  true  DOA. 

Figure  1 1 .2  shows  the  Root  Mean  Square  Error  (RMSE)  of  the  DOA  estimation  of  a 
signal  arriving  from  45°  using  the  ML  and  IMUSIC  methods  for  different  SNRs.  In  this 
example,  an  8  sensor  ULA  with  element  spacing  A/2  is  used.  Two  uncorrelated  signals  are 
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incident  from  35°  and  45°  and  100  snapshots  are  collected.  The  SNR  is  varied  from  5dB  to 
30dB  and  30  Monte  Carlo  simulations  are  made  for  each  simulation.  From  Figure  1 1.2  it  is 
clear  that  the  RMSE  of  the  IMUSIC  method  is  higher  than  that  of  the  ML  method  but 
approaches  it  asymptotically  for  high  SNRs. 

5 

4.5 


f  3 
!« 

1, 


SNR(dB) 

Fig. 1 1.2  The  RMSE  of  the  DOA  estimation  of  0  =  45°  .for  ML  and  IMUSIC  methods  versus  SNR  using  a  ULA  of 
8  sensors  with  signals  incident  from  35°  and  45°  taking  100  snapshots.  30  Monte  Carlo  simulations  are  made. 


11.2  DOA  and  Mutual  Coupling  Estimation  with  Unknown  Frequencies. 

For  the  DOA  and  mutual  coupling  estimation  with  multiple  unknown  frequencies  the 
Separable  Sub-space  Method  described  in  section  10.2.3  was  used.  A  UCA  with  8  sensors  is 
used  for  the  simulation  and  four  equal-power  narrow-band  sources  centered  with  two  different 
signal  frequencies  are  located  at  the  far  field.  The  four  signals  have  the  following  frequencies 
and  DO  As: 

S(0n,  ft)  =  (30°,  8.5  MHz),  S(0i2  ,fi)  =  (60°,  8.5  MHz), 

S(02i,  f2)  =  (90°,  6.5  MHz),  S(022,  f2)  =  (120°,  6.5  MHz) 

Additive  noise  is  injected  with  a  SNR  of  15  dB  and  100  snapshots  of  array  data  are 
accumulated. 

Figure  1 1.3  shows  the  frequency  estimation.  For  each  of  the  2  estimated  frequencies, 
Figure  1 1.4  and  Figure  1 1.5  show  the  spatial  spectrums  of  the  DOA  estimations  with  both 
unknown  and  estimated  MCMs. 
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Fig. 1 1.3  Frequency  estimation  by  separable  dimension  method  using  a  UCA  of  8  sensors.  Four  signals  with 
frequencies  (8.5  MHz,  6.5  MHz)  are  incident  from  (30°,  60°)  and  (90°,  120°)  with  100  snapshots.  The  SNR  =15  dB. 


drecbon  cA  arrival  (degree) 

Fig.  1 1 .4  The  spatial  spectrum  for  estimated  MCM  and  unknown  MCM  at  the  estimated  frequency  fj  =  8.5  MHz 
using  a  UCA  of  8  sensors.  Four  signals  with  frequencies  (8.5  MHz,  6.5  MHz)  are  incident  from  (30°,  60°)  and  (90°,  120°) 
with  100  snapshots.  The  SNR  =15  dB. 
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Fig.  1 1 .5  The  spatial  spectrum  for  estimated  MCM  and  unknown  MCM  at  the  estimated  frequency  f2  =  6.5  MHz 
using  a  UCA  of  8  sensors.  Four  signals  with  frequencies  (8.5  MHz,  6.5  MHz)  are  incident  from  (30°,  60°)  and  (90°,  120°) 
with  100  snapshots.  The  SNR  =15  dB. 

By  inspecting  the  simulation  results  in  Figures  1 1 .4  and  1 1 .5,  it  can  be  seen  that  the 
mutual  coupling  effect  does  not  greatly  affect  the  frequency  estimation  but  severely  affects  the 
DOA  estimation  when  using  the  separable  dimension  estimation  method. 


Fig.  1 1.6  The  contours  of  the  2D  MUSIC  spectrum  for  DOA  and  frequency  with  compensated  MCM  where  the 
dotted  lines  represent  the  true  DOA  and  frequency  pairs. 


Compensated  with  the  estimated  MCM,  the  separable  dimension  estimation  method  can 
estimate  DOAs  that  are  very  close  to  the  true  DOAs.  Figure  11.6  shows  the  contours  of  the 
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2D  MUSIC  spectrum  for  DO  A  and  frequency  with  compensated  MCM.  The  dotted  lines 
indicate  the  true  DOA/frequency  pairs. 


Fig.  1 1.7  The  MSE  of  DO  A  estimation  of  9  =30°  using  separable  dimension  method  versus  SNR  using  a  UCA  of  8 
sensors  and  30  Monte  Carlo  simulations.  Four  signals  with  frequencies  (8.5  MHz,  6.5  MHz)  are  incident  from  (30°,  60°)  and 
(90°,  120°)  with  100  snapshots. 


Fig.  1 1.8  The  MSE  of  DO  A  estimation  of  f  =  8.5  MHz  using  separable  dimension  method  versus  SNR 
using  a  UCA  of  8  sensors  and  30  Monte  Carlo  simulations.  Four  signals  with  frequencies  (8.5  MHz,  6.5  MHz)  are  incident 
from  (30°,  60°)  and  (90°,  120°)  with  100  snapshots. 


Figures  1 1.7  and  1 1.8  illustrate  the  mean  square  error  (MSE)  of  the  estimated  DOA  and 
frequency,  DOA  =30°  and  f  =  8.5  MHz,  and  also  compare  them  to  their  theoretical  lower 
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bound.  In  this  example,  the  array  and  source  scenarios  are  the  same  as  the  above  example 
except  that  150  snapshots  are  taken.  The  SNR  is  varied  from  5  dB  to  35  dB  and  30  Monte 
Carlo  simulations  are  made  for  each  test. 


11.3  Array  Sensor  Gain/Phase  Error  Calibration. 

To  simulate  the  simultaneous  estimation  of  frequency,  DO  A  and  sensor  gain/phase 
error  the  Alternating  Iterative  method  of  section  10.2.4  and  the  Global  Optimisation  using 
Simulated  Annealing  of  section  10.2.5  were  considered.  For  both  simulations  a  UCA  of  8 
omnidirectional  sensors  is  used. 

11.3.1  Alternating  Iterative  Method 

In  this  example  four  equal-power  narrow-band  sources  are  located  in  the  far  field  of  the 
array  with  the  frequencies  and  directions: 

S(0u,  fi)  =  (30°,  7.5  MHz),  S(012  ,fi)  =  (58°,  7.5  MHz), 

S(02i,  f2)  =  (85°,  5.5  MHz),  S(022,  f2)  =  (1 15°,  5.5  MHz) 

The  SNR  was  set  at  20  dB  and  100  data  samples  were  accumulated. 


DOA  esdamt«3n(apatlal  spectrun)  for  f  *  7.5  MHz  after  iterations  1.  5. 10 


Fig. 11.9  MUSIC  spectrum  for  the  Alternating  Iterative  method  after  iterations  1,6, and  10  using  a  UCA  of  8 
sensors  with  two  signals  with  frequency  f  i  =7  .5  MHz  incident  from  30°  and  58°.  100  snapshots  with  SNR  =20  dB 
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Figures  1 1.9  and  11.10  show  the  MUSIC  spectrums  after  iterations  1,  6  and  10  for  the 
different  estimated  frequencies,  f  =  [7.5,  5.5]  MHz,  respectively.  It  is  clear  from  the  results 
that  after  10  iterations,  the  array  sensor  gain/phase  errors  have  been  corrected  and  the  estimated 
DOAs  are  close  to  the  true  DOAs. 


DQAesUamfcon  (spatial  spectrum )  lor  ( •>  5.5  MHz  after  iterations  1,5,10 


Fig. 11. 10  MUSIC  spectrum  for  the  Alternating  Iterative  method  after  iterations  1,6, and  10  using  a  UCA  of  8 
sensors  with  two  signals  with  frequency  f  2  =5  .5  MHz  incident  from  85°  and  115°.  100  snapshots  with  SNR  =20  dB 


11.3.2  Simulated  Annealing  Method. 

Using  the  same  array  structure  as  the  previous  example,  this  uses  a  Simulated 
Annealing  algorithm  to  estimate  the  frequency  and  DOA  of  four  signals  in  the  presence  of 
sensor  gain/phase  errors.  Four  equal-power  narrow-band  sources  are  located  in  the  far  field  of 
the  array  with  the  frequencies  and  directions: 

S(0n,  fi)  =  (15°,  5.5  MHz),  S(0i2  ,fi)  =  (120°,  5.5  MHz), 

S(02i,  f2)  =  (150°,  7.5  MHz),  S(022,  f2)  =  (250°,  7.5  MHz) 

The  SNR  =  10  dB  and  100  samples  are  accumulated.  Since  the  annealing  schedule 
(To,  no)  ...  (Tf,  nf)  is  a  critically  important  element  in  the  simulated  annealing  algorithm,  the 
annealing  schedule  should  be  selected  carefully  to  guarantee  that  the  cost  function  converges  to 
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the  global  optimum.  Note  that  for  different  optimization  problems,  the  annealing  schedule 
selection  will  be  different.  The  annealing  schedule  consists  of: 

Initial  temperature  To:  A  value  To  should  is  chosen  that  enables  the  system  to  climb 
out  of  any  local  minimum. 

Equilibrium  condition  n  j :  The  equilibrium  condition  is  met  when  a  certain  number  of 
inner  circulations  are  accepted  at  that  temperature. 

Temperature  decrease  T;+i:  When  each  equilibrium  state  is  reached,  the  temperature 
will  be  decreased  by  a  constant  value. 

Stop  criteria  Tp  A  certain  number  of  outer  circulations.  Also,  a  pre-determined 
threshold  may  select  it;  when  the  cost  function  at  the  end  of  2  consecutive  accepted  values  is 
less  than  this  number,  the  outer  circulation  will  stop  and  the  cost  function  is  assumed  to  have 
reached  its  minimum. 

Figures  11.11  and  11.12  show  the  results  of  the  DOA  and  frequency  estimations  using 
the  simulated  annealing  method.  As  can  be  seen,  with  random  searching  in  DOA  and 
frequency  space  by  the  simulated  annealing  algorithm,  the  curves  finally  converge  to  the 
corresponding  DOAs  and  frequencies  and  the  DOAs  and  frequencies  are  estimated 
simultaneously.  8000  annealing  iterations  are  shown  in  the  graphs.  Figure  11.13  shows  the 
cost  function  for  the  simulated  annealing  method  where  the  cost  function  curve  moves  to  the 
minimum  from  local  stationary  points.  This  allows  the  algorithm  to  proceed  and  to  find  a 
globally  optimal  solution. 
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Fig.  1 1.1 1  DOA  estimation  using  simulated  annealing  A  UCA  of  8  sensors  is  used.  4  signals  with  frequencies,  f[  = 
7.5  MHz  and  f2  =  5.5  MHz  are  incident  from  (150°,  250°)  and  (15°,  120°).  100  snapshots  with  SNR  =10dB. 
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Fig. 11. 12  Frequency  estimation  using  simulated  annealing  A  UCA  of  8  sensors  is  used.  4  signals  with 
frequencies,  fi  =  7.5  MHz  and  f2  =  5.5  MHz  are  incident  from  (150°,  250°)  and  (15°,  120°).  100  snapshots  with  SNR  =10dB. 


The  number  of  iterations 


Fig.  1 1 . 1 3  The  cost  function  for  the  simulated  annealing  method 
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12.  Conclusions. 

12.1  Introduction. 

This  report  is  the  final  report  of  a  project  undertaken  by  the  University  of  Limerick  on 
behalf  of  the  Air  Force  Research  Laboratory  (AFRL/IFGA,  Rome,  NY)  through  the  European 
Office  of  Aerospace  Research  and  Development.  The  main  objective  of  the  project  is  to 
investigate  the  possibility  of  using  the  known  characteristics  of  GPS  signals  to  calibrate  smart 
antenna  arrays  ‘on-the-fly’  and  thus  reduce  the  time  taken  to  acquire  and  characterise  other 
unknown  signals,  especially  in  hostile  environments.  This  necessitates  the  development  of 
improved  robust,  stable  and  accurate  estimation  algorithms  to  process  the  data  sampled  from 
the  sensors  of  an  antenna  array.  In  practice,  it  may  be  advantageous  to  utilise  different 
algorithms  to  handle  different  parts  of  the  problem. 

12.2  Antenna  array  model  and  systemic  errors. 

In  section  7  it  was  shown  that,  for  an  ideal  antenna  array,  the  azimuth  and  elevation  of 
an  unknown  signal  of  known  frequency  can  be  calculated  from  the  measured  time  delays  at 
each  sensor  by  reference  to  the  ideal  antenna  manifold;  which  is  itself  a  function  of  the 
individual  sensor  response  and  the  array  geometry.  In  practice,  however,  the  array  manifold  is 
distorted  by  gain  and  phase  errors  between  the  sensors  and  their  associated  circuitry  and  mutual 
coupling  between  array  elements.  The  gain  and  phase  errors  of  the  sensors  are  measured  with 
respect  to  a  reference  sensor  and  are  independent  of  each  other;  they  can  be  described  as  a 
diagonal  matrix.  The  effects  of  mutual  coupling  are  more  complex  as  they  are  dependent  on 
the  position  of  the  sensor  within  the  array.  The  effects  are  best  described  using  a  matrix,  the 
Mutual  Coupling  Matrix  (MCM),  which  in  the  usual  case  of  symmetric  arrays  has  a  symmetric 
structure;  for  a  ULA  the  MCM  is  a  banded  Toeplitz  matrix  and  for  a  UCA  it  is  a  circulant 
matrix. 

In  particular,  these  errors  are  frequency  dependent  and  thus  severely  degrade  the 
accuracy  of  spatial  or  spatio-temporal  parameter  estimations  when  handling  multiple  signals 
with  different  frequencies.  In  this  study,  therefore,  much  effort  has  been  put  into  investigating 
and  proposing  methods  of  estimating  the  gain/phase  and  mutual  coupling  errors  in  the  situation 
of  both  known  and  unknown  signal  frequencies. 
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12.3  Parameter  estimation  methods. 

For  the  case  of  signals  with  known  frequency  (as  would  be  the  case  when  using  GPS 
signals)  a  Maximum  Likelihood  method  (section  10.2.1)  and  an  Iterative  MUSIC  method 
(which  is  a  Subspace  method)  (section  10.2.2)  were  considered  for  the  estimation  of  directions 
of  arrival  (DOAs)  and  the  MCM  simultaneously.  Analysis  of  the  effects  of  model  errors  and 
sensitivity  analysis  (section  9)  appeared  to  show  that  the  Subspace  method  would  be  more 
accurate.  However,  when  simulations  were  run,  the  ML  method  showed  a  faster  convergence 
to  the  actual  value  (section  11.1)  and  a  lower  RMSE. 

For  the  case  of  multiple  unknown  frequency  signals  a  Separable  Dimension  Subspace- 
based  method  was  considered  (section  10.2.3)  to  simultaneously  estimate  the  signal  frequencies 
and  DOAs  as  well  as  the  mutual  coupling  parameters.  With  separable  dimension  processing 
the  spatial  and  temporal  estimation  problems  are  separated;  the  signal  frequencies  are  first 
estimated  using  a  subspace  method  and  then  the  DOAs  and  MCM  are  estimated  at  each  of  the 
estimated  frequencies.  In  this  way  the  computational  complexity  is  relatively  low  and  the 
perfonnance  is  satisfactory.  Two  separable  subspace  algorithms  were  developed  to 
demonstrate  the  power  of  this  approach  -  an  Alternating  Iterative  method  and  a  Global 
Optimisation  method  based  on  Simulated  Annealing.  This  estimation  procedure  is  a  non¬ 
linear,  multi-dimensional  optimisation  problem  and  whether  or  not  the  algorithm  can  converge 
to  a  global  optimum  depends  on  the  choice  of  initial  values.  When  simulations  were  run 
(section  1 1.3),  both  algorithms  were  shown  to  be  capable  of  accurate  frequency  and  DOA 
estimation  with  effective  correction  of  systematic  errors.  The  Global  Optimisation  method  is 
far  more  computationally  intensive  as  well  as  being  sensitive  to  the  initial  annealing 
parameters;  however,  it  does  have  the  advantage  of  ensuring  a  high  probability  of  convergence 
to  a  global  optimum  whereas  this  cannot  be  guaranteed  using  the  Alternating  Iterative  method. 

12.4  The  use  of  GPS  signals  to  estimate  the  MCM  and  gain/phase  errors. 

Signals  from  GPS  satellites  are  ideal  for  calibrating  an  antenna  array  because  they  are  of 
the  same  known  frequency  and,  once  the  position  and  attitude  of  the  array  is  known,  the  DOA 
of  each  signal  is  known  accurately  and  independent  of  array  parameters.  This  would  enable  an 
ML  algorithm,  such  as  that  developed  in  section  10.2.1,  to  rapidly  estimate  the  MCM  and 
gain/phase  errors;  either  as  separate  matrices  or  as  an  actual  array  manifold.  It  might  be 
advantageous  to  carry  out  this  calibration  at  both  LI  and  L2  so  as  to  obtain  some  information 
on  the  frequency  dependence  of  the  error  matrices. 
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Once  we  have  obtained  accurate  MCM  and  gain/phase  error  matrices,  they  may  be  used, 
with  rough  frequency  correction  if  available,  as  initial  conditions  for  using  a  separable 
dimension  algorithm  (such  as  derived  in  sections  10.2.4  or  10.2.5)  to  estimate  the  frequency 
and  DOA  of  unknown  signals.  Simulations  in  section  1 1  have  shown  that  the  speed  of 
convergence  and  accuracy  of  these  estimation  methods  is  very  dependent  on  the  initial 
conditions.  The  use  of  these  improved  initial  conditions  should,  therefore,  speed  up  the 
estimation  of  the  unknown  parameters  and  increase  the  chance  of  attaining  a  global  optimum. 
This  may  enable  the  faster  Alternating  Iterative  method  to  be  used  instead  of  the  alternative 
Global  Optimisation  using  Simulated  Annealing.  Alternatively  it  may  lead  to  faster 
convergence  using  Simulated  Annealing;  further  investigation  is  necessary  to  evaluate  the 
advantages  of  more  accurate  initial  conditions. 

We  have  shown  that  GPS  signals  may  be  used  to  calibrate  an  antenna  array  by 
estimating  the  MCM  and  gain/phase  error  matrix.  This  would  provide  improved  initial 
conditions  for  algorithms  used  to  estimate  the  frequency  and  DOA  of  unknown  signals,  thus 
ensuring  faster  and  more  accurate  estimations. 


U.  L.  Antenna  Research  Group 


Page  63  of  75 _ |  Document  ID:  ARG-ARFL-PJ1.0/  VER  1.7 


Air  Force  Research  Laboratory  -EOARD 


Contract  No.F73001  F30602  99MV072 


13.  Bibliography. 

[1]  H.  Krim  and  M.  Viberg,  Two  Decades  of  Array  Signal  Processing  Research:  The  Parametric 
Approach,  IEEE  Signal  Processing  Magazine,  13(4):  67  -  94,  July  1996. 

[2]  G.  V.  Tsoulos,  Smart  Antennas  for  Mobile  Communications  Systems;  benefits  and 
Challenges,  IEE  Electronics  and  Communication  Engineering  Journal,  84  -  94,  April  1994. 

[3]  M.  S.  Bartlet,  Smoothing  Periodograms  from  Time  Series  with  Continuous  Spectra, 

Nature,  Vol.  161:686  -  687,  1948. 

[4]  J.  Capon,  High-resolution  Frequency-wavenumber  Spectrum  Analysis,  Proc.  of  IEEE,  57(8): 
1408  -  1418,  August  1969. 

[5]  R.  O.  Schmidt,  Multiple  Emitter  Location  and  Signal  Parameter  Estimation,  IEEE  Trans,  on 
Antenna  and  Propagation,  34(3):  276  -  280,  March  1986. 

[6]  M.  Viberg  and  B.  Ottersten,  Sensor  Array  Processing  Based  on  Subspace  Fitting,  IEEE 
Trans,  on  Signal  Processing,  39(5):  1110  -  1121,  May  1991. 

[7]  R.  Roy  and  T.  Kaliath,  ESPRIT  -  Estimation  of  Signal  Parameters  via  Rotational  Invariance 
Techniques,  IEEE  Trans,  on  Acoustic,  Speech  and  Signal  Processing,  37(7):984  -  995,  1989. 

[8]  R.  A.  Fisher,  On  an  Absolute  Criterion  for  Fitting  Frequency  Curves,  Messenger  of 
Mathematics,  Vol.  41:  155  -  160,  1912. 

[9]  B.  Friedlander  and  A.  J.  Weiss,  Direction  Finding  in  the  Presence  of  Mutual  Coupling,  IEEE 
Trans,  on  Antennas  and  Propagation,  39(3):  273  -284,  March  1991 

[10]  J.  Mao,  B.  Champagne,  M.  O’Droma  and  K.  Kwiat,  Separable  Dimension  Subspace  Method 
for  the  Joint  Signal  Frequencies,  DOAs  and  Sensor  Mutual  Coupling  Estimation, 

[11]  H.  Akaike,  Fitting  Autoregressive  Models  for  Prediction,  Ann.  Inst.  Statist.  Math.,  Vol. 
21,  pp.243  -247,  1969 

[12]  P.  Stoica  and  R.  Moses,  Introduction  to  Spectral  Analysis,  Prentice  Hall,  1997 

[13]  E.  Aartts  and  J.  Korst,  Simulated  Annealing  and  Boltzmann  Machines,  John  Wiley  &  Sons, 
1990. 

[14]  S.  Geman  and  G.  Geman,  Stochastic  Relaxation,  Gibbs  Distribution,  and  Bayesian 
Restoration  of  Image,  IEEE  Trans,  on  PAMI,  Vol.6:  721  -  741,  1984. 

[15]  R.  Ahmed  and  B.  L.  Evans,  Optimization  of  Signal  Processing  Algorithms,  Proc.  of  IEEE  , 
pp. 1401  -  1406,  1997. 

[16]  C.  A.  Balanis,  Antenna  Theory:  Analysis  and  Design,  John  Wiley  &  Sons,  Inc.,  NY  1982. 


U.  L.  Antenna  Research  Group 


Page  64  of  75 _ |  Document  ID:  ARG-ARFL-PJ 1 .0/  VER  1.7 


Air  Force  Research  Laboratory  -EOARD 


Contract  No.F73001  F30602  99MV072 


[17]  G.  C.  Brown,  J.  H.  McClellan  and  E.  J.  Holder,  Eigenstructure  Approach  for  Array 
processing  and  Calibration  with  General  Phase  and  Gain  Perturbations,  Proc.  of  IEEE 
ICASSP'91 ,  pp.  1365  -  1368,  1991. 

[18]  S.  Chen,  R.  H.  Istepanian,  and  B.  L.  Luk,  Signal  Processing  Applications  Using  Adaptive 
Simulated  Annealing,  Proc.  of  IEEE  ,  pp.  842  -  849,  1999. 

[19]  Y.  Chen  and  C.  Chen,  Direction-of-arrival  and  Frequency  Estimations  for  Narrowband 
Sources  Using  Two  Single  Rotation  Invariance  Algorithms  with  the  Marked  Subspace,  Radar 
and  Signal  Processing,  IEE  Proceedings  -  F,  139  (4)  :  297  -300,  August  1992. 

[20]  M.  P.  Clark,  and  L.  L.  Scharf,  Two-dimensional  Modal  Analysis  Based  on  Maximum 
Likelihood,  IEEE  Trans,  on  Signal  Processing,  42(6):  1443  -  1452,  June  1994. 

[21]  H.  Cramer,  Mathematical  Methods  of  Statistics,  Princeton  University  Press,  1946. 

[22]  A.  D.  Craig,  C.  K.  Leong  and  A.  T.  Wishar,  Digital  Signal  Processing  in  Communications 
Satellite  Payloads,  Electronics  and  Communications  Engineering  Journal,  pp.107  -  114,  June 
1992. 

[23]  R.  S.  Elliot,  Antenna  Theory  and  Design,  Prentice-Hall,  Englewood  Cliffs,  NJ,  1983. 

[24]  Emile  Aartts  and  Jan  Korst,  Simulated  Annealing  and  Boltzmann  Machines,  John  Wiley  \& 
Sons,  1990. 

[25]  R.  B.  Ertel,  P.  Cardieri,  K.  W.  Sowerby,  T.  S.  Rappaport  and  J.  H.  Reed,  Overview  of 
Spatial  Channel  Models  for  Antenna  Array  Communication  Systems,  IEEE  Personal 
Communications,  pp.10  -  22,  February  1998. 

[26]  B.  Friedlander,  A  Sensitivity  Analysis  of  the  MUSIC  Algorithm,  IEEE  Trans,  on 
Acoustics,  Speech,  and  Signal  Processing,  38(10),  October  1990. 

[27]  A.  J.  Weiss  and  B.  Friedlander,  Eigenstrecture  Methods  for  Direction  Finding  with  Sensor 
Gain  and  Phase  Uncertainties,  Circuits,  Systems  and  Signal  Processing,  Vol.  9,  No.  3,  pp.272  - 
300,  1990. 

[28]  E.  M.  Friel  and  K.  M.  Pasala,  Wideband  Bearing  Estimation  with  Compensation  for 
Mutual  Coupling  Effects,  Proc.  of  IEEE  Antenna  and  Propagation  Symposium  (AP-S/URST94)  , 
Vol.  3,  pp.1556-  1559,  1994. 

[29]  S.  Geman  and  G.  Geman,  Stochastic  Relaxation,  Gibbs  Distribution,  and  Bayesian 
Restoration  of  Image,  IEEE  Trans,  on  PAMI,  Vol.6:  721  -  741,  1984. 

[30]  I.  J.  Gupta  and  A.  K.  Ksienski,  Effect  of  Mutual  Coupling  on  the  Perfonnance  of  Adaptive 
Arrays,  IEEE  Trans,  on  Antennas  and  Propagation,  31(5):  785-791,  September  1983. 

[31]  M.  H.  Hayes,  Statistical  Digital  Signal  Processing  and  Modeling,  John  Wiley  and  Sons, 
Inc. 


Air  Force  Research  Laboratory  -EOARD 


Contract  No.F73001  F30602  99MV072 


[32]  S.  Haykin,  Advances  in  Spectrum  Analysis  and  Array  Processing,  Vol.  2,  Prentice-Hall, 
1991. 

[33]  M.  Hamid  and  R.  Mamid,  Equivalent  Circuit  of  Dipole  Antenna  of  Arbitrary  Length, 
IEEE  Trans,  on  Antennas  and  Propagation,  45(11):  1695  -  1696,  November  1997. 

[34]  A.  S.  Householder,  The  Theory  of  Matrices  in  Numerical  Analysis,  NEW  York,  1975. 

[35]  D.  H.  Johnson  and  D.  E.  Dudgen,  Array  Signal  Processing-Concepts  and  Techniques, 
Prentice-Hall,  1993. 

[36]  S.  Kirkpatric,  C.  D.  Gellat  and  M.  Vecchi,  Optimization  by  Simulated  Annealing, 
Science, V ol.  220,  pp.671  -  680,  1983 

[36]  A.  Leshem  and  M.  Wax,  Array  calibration  in  the  presence  of  multipath,  IEEE  Trans,  on 
Signal  Processing,  48(1):  53  -  59,  Jan.  2000. 

[37]  J.  T.  -H.  Lo  and  S.  L.  Marple,  Observability  Conditions  for  Multiple  Signal  Direction 
Finding  and  Array  Sensor  Localization,  IEEE  Trans,  on  Signal  Processing,  Vol. 40:  2641  - 
2650, 1992. 

[38]  S.  Lundgren,  A  Study  of  Mutual  Coupling  Effects  on  the  Direction  Finding  Performance 
of  ESPRIT  with  a  Linear  Microstrip  Patch  Array  Using  the  Method  of  Moments,  IEEE  AP-S 
Symposium,  pp.1372  -  1375,  Baltimore,  USA,  July  1996. 

[39]  D.  G.  Manolakis,  V.  K.  Ingle,  and  S.  M.  Kogon,  Statistical  and  Adaptive  Signal 
Processing,  McGraw-Hill ,  2000. 

[40]  J.  Mao  and  M.  Viberg,  Joint  Estimation  for  Frequencies,  Bearings  and  Array  Model  Errors 
by  Using  Simulated  Annealing  ,  IASTED  Conference  on  Modelling  and  Simulation  (MS'98), 
Pittsburgh,  USA,  May  1998. 

[41]  B.  C.  Ng  and  C.  M.  S.  See,  Sensor-Array  Calibration  Using  a  Maximum-Likelihood 
Approach  ,  IEEE  Trans,  on  Antennas  and  Propagation,  44(6):  827  -  835,  June  1996. 

[42]  B.  Ottersten,  M.  Viberg,  P.  Stoica,  and  A.  Nehorai,  Exact  and  large  Sample  ML 
Techniques  for  Parameter  Estimation  and  Detection  in  Array  Processing,  Radar  Array 
Processing },  pp.99  -  151,  (edited  by  Haykin,  Litiva,  and  Shepherd),  Springer- Verlag,  Berlin, 
1993. 

[43]  K.  M.  Pasala  and  E.  M.  Friel,  Mutual  Coupling  Effects  and  Their  Reduction  in  Wideband 
Direction  of  Arrival  Estimation,  IEEE  Trans,  on  Aerospace  and  Electronic  Systems),  30(4): 
1116-  1122,  October  1994. 

[44]  A.  J.  Pauraj  and  T.  Kailath,  Direction  of  Arrival  Estimation  by  Eigenstructure  Methods 
with  Unknown  Sensor  Gain  and  Phase,  Proc.  ofICASSP-85 :  17.7.1  -  17.7.3,  1985. 


U.  L.  Antenna  Research  Group 


Page  66  of  75 _ |  Document  ID:  ARG-ARFL-PJ 1 .0/  VER  1.7 


Air  Force  Research  Laboratory  -EOARD 


Contract  No.F73001  F30602  99MV072 


[45]  A.  J.  Paulraj,  R.  O.  Roy,  and  T.  Kaliath,  A  Subspace  Rotation  Approach  to  Signal 
Parameter  Estimation,  Proc.  of  IEEE,  74(7):  1044-  1045,  1986. 

[46]  D.  T.  Pham  and  D.  Karaboga,  Intelligent  Optimization  Techniques,  Springer- Verlag, 
2000. 

[47]  J.  Pierre  and  M.  Kaveh,  Experimental  Performance  of  Calibration  and  Direction-Finding 
Algorithms,  IEEE  ICASSP'91 ,  Vol.2,  pp.1365  -  1368,  Toronto,  Canada,  May  1991. 

[48]  B.  A.  Pontano,  Satellite  Communications:  Services,  Systems  and  Technologies,  IEEE 
MTT-S  Digest,  pp.l  -  4,  1998. 

[49]  K.  C.  Sharman,  Maximum  Likelihood  Parameter  Estimation  by  Simulated  Annealing, 
Proc.  of  ICASSP-88,  pp.274 1  -  2744,  1988. 

[50]  I.  S.  D.  Solomon,  Y.  I.  Abramovich,  D.  A.  Gray,  and  S.  J.  Anderson,  Perfonnance  of 
OTH  Radar  Array  Calibration,  IEEE  ICASSP'98,  Vol.4,  pp.2025  -  2028,  Seattle,  USA,  May 
1998. 

[51]  I.  S.  D.  Solomon,  D.  A.  Gray,  Y.  I.  Abramovich,  and  S.  J.  Anderson,  Over-the-horizon 
Radar  Array  Calibration  Using  Echoes  from  Ionized  Meteor  Trails,  IEE  Proc.  -Radar,  Sonar 
Navig.,  Vol.145,  No. 3,  pp.173  -  180,  June  1998. 

[52]  Peter  Stoica  and  Arye  Nehorai,  MUSIC,  Maximum  Likelihood,  and  Cramer-Rao  Bound, 
IEEE  Trans,  on  Acoustics,  Speech  and  Signal  Processing,  37(5):  720  -741,  May  1989. 

[53]  T.  Svantesson,  The  Effects  of  Mutual  Coupling  on  the  Direction  Finding  Accuracy  of  a 
Linear  Array  of  Thin  Dipoles  of  Finite  Length,  Technical  report  R007/1998,  Department  of 
Signals  and  Systems,  Chalmers  University  of  Technology,  Goteborg,  Sweden,  August  1998. 

[54]  A.  Swindlehurst  and  T.  Kaliath,  A  Performance  Analysis  of  Subspace-based  Methods  in 
the  Presence  of  Model  Errors  -  Part  I:  The  Music  Algorithm,  IEEE  Trans,  on  Signal  Processing, 
40(7):  1758  -  1774,  July  1992. 

[54]  A.  Swindlehurst  and  T.  Kaliath,  A  Perfonnance  Analysis  of  Subspace-based  Methods  in 
the  Presence  of  Model  Errors  -  Part  II:  Multidimensional  Algorithms,  IEEE  Trans,  on  Signal 
Processing,  41(9):  2882  -2890,  September  1993. 

[55]  M.  Viberg  and  A.  L.  Swindlehurst,  A  Bayesian  Approach  to  Auto-Calibration  for 
Parametric  Array  Signal  Processing,  IEEE  Trans,  on  Signal  Processing,  42(12):  3495  -  3506, 
December  1994. 

[56]  M.  Wax  and  I.  Ziskind,  Detection  of  the  Number  of  Coherent  Signals  by  the  MDL 
Principle  ,  IEEE  Trans,  on  Acoustics,  Speech  and  Signal  Processing,  37(8):  pp.l  190  -  1196, 
Aug.  1989. 


U.  L.  Antenna  Research  Group 


Page  67  of  75 _ |  Document  ID:  ARG-ARFL-PJ 1 .0/  VER  1.7 


Air  Force  Research  Laboratory  -EOARD 


Contract  No.F73001  F30602  99MV072 


[57]  S.  R.  White,  Concepts  of  Scale  in  Simulated  Annealing,  Proc.  of  IEEE  Conf.  Computer 
Design:  VLSL  in  Computers,  ICCD-84 ,  pp.646  -  651,  1984 

[58]  H.  Yuan,  K.  Hirasawa,  and  Y.  Zhang,  The  Mutual  Coupling  and  Diffraction  effects  on  the 
Performance  of  a  CMA  Adaptive  Array,  IEEE  Vehicular  Technology,  47(3):  728  -  736,  August 
1998. 

[59]  M.  Zhang  and  Z.  Zhu,  A  Method  for  direction  Finding  under  Sensor  Gain  and  Phase 
Uncertainties,  IEEE  Trans,  on  Antennas  and  Propagation,  vol.  43,  pp.  880  -  883  ,  August  1995. 

[60]  I.  Ziskind  and  M.  Wax,  Maximum  Likelihood  Localization  of  Diversely  Polarized 
Sources  by  Simulated  Annealing,  IEEE  Trans,  on  Antennas  and  Propagation,  38(7):  1111  - 
1114,  July  1990. 


U.  L.  Antenna  Research  Group 


Page  68  of  75 


Document  ID:  ARG-ARFL-PJ1.0/  VER  1.7 


Air  Force  Research  Laboratory  -EOARD 


Contract  No.F73001  F30602  99MV072 


Appendix  A 


Paper  published  in 

Proceedings  of  The  IEEE  34th  Asilomar  Conference  on  Signals,  Systems  and  Computers.  IEEE  o- 
7803-6514-3/00,  PP605-609. 


October  29  -  November  1,  2000,  in  Pacific  Grove,  California. 

J.  Mao,  B.  Champagne,  M.  O’Droma  &  K.  Kwiat.  Separable  dimension  subspace  method  for 
joint  signal  frequencies,  DOAs  and  sensor  mutual  coupling  estimation’. 


U.  L.  Antenna  Research  Group 


Page  69  of  75 _ |  Document  ID:  ARG-ARFL-PJ 1 .0/  VER  1.7 


Air  Force  Research  Laboratory  -EOARD 


Contract  No.F73001  F30602  99MV072 


Proceedings  of 

The  IEEE  34th  Asilomar  Conference  on  Signals,  Systems  and  Computers 

October  29  -  November  1,  2000,  in  Pacific  Grove,  California. 


Separable  Dimension  Subspace  Method  for  Joint  Signal  Frequencies,  DOAs  and 
Sensor  Mutual  Coupling  Estimation 


Jian  Mao1,2,  Benoit  Champagne1,  Mairtin  O’Droma2,  and  Kevin  Kwiat3 


1  Dept,  of  Electrical  and  Computer  Engineering 
McGill  University,  Montreal,  Quebec  H3A  2A7,  Canada 
mao@tsp.ece.mcgill.ca,  champagne@ece.mcgill.ca 

2Dept.of  Electronic  and  Computer  Engineering 
University  of  Limerick,  Limerick,  Ireland  * 
jain.mao,  mairtin.odroma@ul.ie 


3Air  Force  Research  Lab 
AFRL/IFGA,  Rome,  NY  13441-4505,  USA 
kwiatk@rl.af.mil 


Abstract 

To  extract  the  frequencies  and  direction  of  arrivals 
(DOAs)  of  multiple  sources  from  experimental  data  col¬ 
lected  by  a  sensor  array  is  a  multiple  parameter  estimation 
problem.  Some  important  algorithms  for  spatial-temporal 
processing  have  been  developed  in  the  past  decades.  A 
practical  problem,  not  often  considered,  is  that  the  differ¬ 
ent  sensors  in  the  array  affect  each  other  through  mutual 
coupling.  This  effect  varies  with  frequencies  and  degrades 
the  performance  of  algorithms.  Thus,  a  separable  dimen¬ 
sion  subspace  method  to  simultaneously  estimate  signal  fre¬ 
quencies,  direction  of  arrivals  ( DOAs )  and  sensor  mutual 
coupling  is  proposed  in  this  paper. 


1.  Introduction 

The  estimation  of  the  frequencies  and  direction  of  ar¬ 
rivals  (DOA's)  of  multiple  signals  by  using  an  antenna  array 
has  attractive  and  important  applications  in  various  areas, 
such  as  radar  and  communication  systems.  In  particular, 
for  next  generation  mobile  communication  systems,  due  to 
an  endless  quest  for  increased  capacity  and  improved  qual¬ 
ity,  frequency  and  DOA  estimation  become  a  requirement. 
Many  effective  algorithms  for  simultaneously  estimating 

•Thanks  to  American  Air  Force  Research  Lab  for  their  funding  the  re¬ 
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the  frequencies  and  DOA  of  multiple  signals  have  been  de¬ 
veloped  in  the  recent  past  [1,  3, 4],  Most  of  them  assume 
that  the  array  response  is  completely  known.  However, 
many  factors,  such  as  mutual  coupling  among  the  different 
array  sensors,  will  alter  the  array  response  in  practical  ap¬ 
plications.  Particularly,  in  the  presence  of  multiple  sources 
with  different  frequencies,  the  sensor  mutual  coupling  will 
vary  with  frequency  and  thus  severely  affect  frequency  and 
DOA  estimation  accuracy.  In  Reference  [5],  a  joint  signal 
frequencies,  DOAs  and  array  model  errors  estimation  algo¬ 
rithm  was  proposed.  Due  to  simulated  annealing  process 
employed,  the  proposed  algorithm  has  high  computation 
cost.  In  this  paper,  a  separable  dimension  subspace  method 
to  simultaneously  estimate  signal  frequencies,  DOAs  and 
mutual  coupling  parameters  of  antenna  array  is  presented. 
With  separable  dimension  processing,  a  spatial  and  tem¬ 
poral  estimation  problem  is  separated,  i.e,  the  frequencies 
are  first  estimated  by  using  a  subspace  method  and  then  the 
DOAs  and  mutual  coupling  parameters  are  estimated  for  the 
each  estimated  frequency.  In  this  way,  the  computational 
complexity  of  the  proposed  method  is  relatively  small. 

2.  Problem  Formulation 

Consider  an  array  composed  of  M  sensors  and  each  sen¬ 
sor  output  being  fed  to  a  tap-delay-line  with  m  taps.  The 
delay  between  adjacent  taps  is  to-  Specially,  let  X[(t)(i  = 
1, 2,  •  •  • ,  M )  denotes  the  output  of  the  1-th  sensor  at  time  t, 
and  let  X;(t)  =  [xi(f),  •  •  •  ,xj(t  —  (m  —  l)to)]T-  Assume 
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that  p  narrow-band  sources  with  different  frequencies  and 
directions  impinge  on  the  array  and  the  coming  signals  can 
be  divided  into  D  groups  according  to  their  wavelength,  i.e, 

D 

p  =  £pfc  (D  <  p)  (1) 

k= 1 

where  pt  is  the  number  of  signals,  from  different  direc¬ 
tions,  in  the  k  -  th  group,  which  wavelength  is  A*  ( k  — 
1,2, 

The  I-th  sensor  (I  =  1, 2,  -  ■  • ,  M)  output  vector  may  be 
defined  as: 

X,(f) 

D  M 

=  E  +  N,(t); 

fc=li=l 

(fc  =  1,2,  •  •  ■  ,D,i  =  1, 2,  •  ■  • ,  Af) 


3  Subspace  Method 

Assume  that  the  number  of  signals  p  and  the  number  of 
frequency  groups,  D,  are  known. 

Define 

1  M 

A  =  ±  Y,  EMf  -(h~  l)fo)X,H(t)]  (4) 

W  1=1 

(h=  1,2, -••,£>) 

1  m 

tj|  =  -i-  52  E[»lb  -  (/»  -  l)fo)Yw(t  -  (h  -  l)f0)]  (5) 

mtet 

0  =  1,2,  ■••,?) 

where 


where 

bi(0*)  =  [exp  (-Mn(5E,i))---exp(-jtu*ri(5E,pfc))] 
is  a  1  x  pt  spatial  steering  vector; 

“(we)  =  [1  exp  (-jw*to)  •  •  •  exp  (-jtUE(m  -  l)to)]  is 
an  m  x  1  temporal  steering  vector;  Ci,i(uk)  is  the  (I,  i )  -  th 
entry  in  the  mutual  coupling  matrix  (MCM)  which  repre¬ 
sents  the  mutual  coupling  effect  from  other  sensors; 

S*(t)  =  [Sk,i(.t)  ■  ■  ■  SkiPi(t)]T  is  a  pt  x  1  signal  vec¬ 
tor  from  the  k-th  group  of  narrowband  sources  and  N;  = 
[jVj(t)  ■  ■  ■  Ni(t  —  (m  —  l)t0)f  is  an  m  x  1  additive  noise 
vector  for  the  i-th  sensor. 

Therefore,  the  array  data  model  can  be  described  as  fol¬ 
lows: 


y,(t  ~(h-  l)t0)  =  (X,(t))„  =  x,(t  -{h-  l)t0), 

Y (t  -  (h  -  l)to) 

=  [Vi{t  ~(h-  l)*o)jte(t  -  (k  -  l)*o>  ' '  -2/m(*  -  {h  -  l)to)]T 
=  [xr(t  -  {h  -  l)f0)®2(t  -  {h  -  l)«o)  l)t0)]r 

and  E(-)  denotes  statistical  expectation. 

{r/J/?=i  an<l  are  frequency  vector  and  direction 

vector.  Their  range  spaces  spanned  by  the  column  of  these 
vectors  which  are  contained  or  are  equal  to  the  range  spaces 
of  a(to)and  B(0),  respectively.  We  have 

n(r'h)  C  R(a(«)),  R(r j|)  C  R(B(0)) 


X(f) 


[XT,(i)Xr2(f)---XrM(0]r 

D  (3) 

E  (C(wE)B(fli))  ®  a(wE)S*(0  +  N(i) 

k=l 


where 


B(Se)  = 


bi(0t) 
bj k(0k) 


is  a  M  x  Pe  spatial  steering  matrix, 
N{i)=[NT1(()NTI(t)-NT„(()]T  is  Mm  x  1 
additive  white  Gaussian  noise,  Symbol  ®  denotes  Kro- 
necker  product  and  C(we)  (fc  =  1,2,  •  •  • ,  D)  is  M  x  M 
mutual  coupling  matrix.  For  a  uniform  linear  array  or 
circular  array,  C(iue)  is  either  a  Toeplitz  matrix  or  a 
circular  matrix,  respectively  [6]. 

Based  on  the  available  samples  { X ( i ) } ,  the  problem 
is  to  estimate  the  frequencies,  DOAs  and  mutual  coupling 
matrix  simultaneously  with  subspace  methods. 


and 

K(ri)JLR(«(«)),  K(u})J.R(B(S)) 

where  SJ(-)  and  N(-)denote  the  range  space  and  null  space. 
As  is  well  known,  the  subspace  methods  are  based  on  the 
above  geometrical  observations.Therefore,we  may  use  cor¬ 
relation  processing  in  spatial  and  temporal  dimension  re¬ 
spectively  to  get  the  estimates  of  {r'h}?=i  and  {tj[  }f=1  and 
then  compute  their  null  subspaces.  Finally,  the  frequencies, 
DOAs  and  mutual  coupling  matrix  can  be  finally  estimated 
with  subspace  methods  by  searching  corresponding  spaces. 

4  Separable  Dimension  Subspace  Algorithm 

Step  1.  Frequency  Estimation 

(1)  Estimate  of  frequency  vectors  r'h  and  rh  : 

M  S 

n  =  ] (6) 
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(7) 

where 

e*  =  l2_y_°  10- ■•{>], 

(h-i) 

(h  =  1,2,  ••■,£>) 

and  S2 3  is  the  estimate  of  the  noise  variance. 

(2)  Gram-Schmidt  (GS)  orthogonalization  and  formation 
of  temporal  projection  matrix  Pu  : 

From  the  vector  {?/,}£,!,  we  can  get  D  orthogonal 
vectors,  {qit}£=1  via  GS  orthogonalization.  Let  Qw  = 
[qi ,  Q2 ,  ■  -  ,  qo],  then  compute  the  temporal  projection  ma¬ 
trix  Pu  =  I  —  Qi  Qff,  which  spans  the  null  space  of 
{»(«*)}?=!■ 

(3)  Estimate  the  unknown  frequencies  with  the  temporal 
projection  matrix  Pu: 

The  frequencies  ate  estimated  as  the  D  largest 

peaks  of  the  function  P(w)  =  (a^  (iu)Pwa(cj))-1,  search¬ 
ing  over  the  frequency  sector  of  interest. 

Step  2.  Direction  and  Mutual  Coupling  Estimation 

(1) Estimate  of  direction  vectors,  and 

i  m  1  ^ 

V,  =  -E  ^  5Z srr(*-(A-l)*o)Y"(*-(fc-l)to)  (8) 

m  a=i  1  i=i 

Vi  =  {Vi  ~  ?2ej )H,  (9) 

where 

e,  =  [0^010  -0] 

(l-i) 

(i  =  1,2,  ••■,?) 

(2)  Via  Gram-Schmit  orthogonalization  of  {rj;  }f=1 ,  the  p 
orthogonal  vectors,  {Ci}f=i  and  spatial  orthogonal  projec¬ 
tion  matrix  Pa  =  X  -  QaQf  are  obtained,  where  Qs=[Cj 

C.-CJ 

(3)  For  each  frequency  estimated  in  Step  1,  the 

directions  9k  and  mutual  coupling  matrix  Cjaijt)  can  be  es¬ 
timated  with  following  equations  iteratively  : 

(BH{0)CH(uk)PeC(uk)B(9))-1 

C(wi b)  =  (G(wjt)-1w)(w"G(wk)_1w)_l 

where  G(u>*)  =  £  B(0fc)HPaB(0,t)  and  w  = 

[i,o,-,  o]L 

5  Computer  Simulations 

A  uniform  circular  array  with  8  sensors  is  used  for  the 
simulation.  Four  equal-power  narrow-band  sources  are 
located  at  the  far  field  of  the  array  with  the  center  frequen¬ 
cies  and  directions:  S(9u,fi)  —  (30°,8.5MHz), 


S(flu,/i)  =  (60°,8.5MHa),  S(S21,/2)  = 

(90°,  6.5 MHz),  S(fi 22,  A)  =  (120°,6.5MiTz).  Ad¬ 
ditive  noise  is  injected  with  SNR  of  15  dB  referenced  to  the 
signal  source.  100  snapshots  of  array  data  are  accumulated. 
Fig.l  shows  the  frequency  estimation.  For  the  2  estimated 
frequencies,  Fig.2  and  Fig.3  show  the  spatial  spectrum 
for  the  DOA  estimations  with  unknown  and  estimated 
MCM,  respectively.  Fig.4  and  Fig5  show  the  mean  square 
error  (MSE)  of  estimated  DOA  and  frequency  and  their 
theoretical  low  bound.  30  Monte  Carlo  simulations  are 
made  for  the  DOA  =  30°  and  /  =  8.5 MHz 


Figure  2.  Spatial  spectrum  at  the  estimated 
frequency  /  =  8 ,5MHz. 


6  Conclusions 

In  this  paper  a  new  algorithm  based  on  separable  dimen¬ 
sion  subspace  method  is  proposed  for  joint  estimation  of 
signal  frequencies,  DOAs  and  the  mutual  coupling  param¬ 
eters  of  antenna  array.  The  presented  algorithm  has  been 
test  by  computer  simulation  studies  and  has  been  found  to 
perform  satisfactorily. 


607 


U.  L.  Antenna  Research  Group 


Page  72  of  75 


Document  ID:  ARG-ARFL-PJ1.0/  VER  1.7 


Air  Force  Research  Laboratory  -EOARD 


Contract  No.F73001  F30602  99MV072 


PROCEEDINGS  OF 

THE  IEEE  34™  ASILOMAR  CONFERENCE  ON  SIGNALS,  SYSTEMS  AND  COMPUTERS 
OCTOBER  29  -  NOVEMBER  1.  2000,  IN  PACIFIC  GROVE,  CALIFORNIA. 


Figure  3.  Spatial  spectrum  at  the  estimated 
frequency  /  =  6.5 MHz. 


Figure  4.  MSE  and  CRLB  for  the  DOA 
estimation,  DO  A  =  30° 


7  Appendix:  Cramer  Rao  Low  Bound 

In  this  appendix,  we  derived  Cramer  Rao  Low  Bound  for 
the  proposed  separable  subspace  method. 

The  likelihood  function  for  the  data  set  is  given  by 

£(X(1),X(2),  —  ,X(JY)/*) 

=  n^rexpt-xo^R-'xy)) 

j- 1  "  ' 

where  X(j)  =  £  C(tu*)B(0t)  0  a(wk)Sk(j)  +  N(j), 
<I>  =  (0,  cu,  C),  and  matrix  R  is 
R  =  E  [XXW] 

=  Y,  C(wfc)B(0t)  ®  a(uk)R,k  (11) 

fc=l 

(C(w*)B(0*)®a(u»*))w+  I 


Figure  5.  MSE  and  CRLB  for  the  frequency 
estimation,  /  =  8.5 MHz 

Assuming  all  sources  have  signal-noise  rate 
P  =  Rsfc /trj,  then  the  Cramer  Rao  low  bound  is 

CRLB($)  =  [J™]-1  (12) 

where  the  m,  n  —  th  element  of  the  Fisher  information  ma¬ 
trix  are 

Jmn  =  N  ■  tr{R_1  ■  3R/d$m  •  R"1  -  3R/a$n}  (13) 

7.1  DOA  terms 

Jen  =  2 N  ■  9t(tr(AePA/fR~1APAefR~I) 

(14) 

+tT(AePAHR-lA9PAHR-1)) 

where  A  =  C  B(0)  ®  a (w)  and  A#  =  C  (9B {6)/d9)  ® 
a(uj)  =  C  B(0)  ®  a(o>). 

7.2  Frequency  to  terms 

Juu  =  2N  ■  $t(tr(AwP AHR~l AP A^R~l) 

(15) 

+fr(Au,PAHR-1A„PAHR-1)) 

where 

Aw  =  C{dB(9)/[du>)  ®  a(w)  +  CB(0)  ®  {da.(u)/dw) 
=  CBu®a  +  CB®£c, 

7.3  Mutual  coupling  terms 

Suppose  circulant  mutual  coupling  matrix  with  only  a 
single  coupling  coefficient  given  by 

Cm  =  ite*  (16) 
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where  p.  and  f  are  magnitude  and  phase  of  mutual  coupling 
coefficient.  Then,  we  have 


=  2N  •  K  (tr  (AflPAHR~lAPA^R~1) 

+  tr(A(1PAwR-1AwFAwR-1)) 

J((  =  2N-X  (tr  (A(PAfiR-lAPAfR-1)  1  ’ 

+  tr  (A^PAWR_1  A{PAhR-1)) 

where  Aa  =  (3C/(9/i)-B<g>a  and  A(  =  (9C/9£)-B<g>a 

7.4  Cross  terms 

The  cross  terms  are  derived  as 

J)u  =  2N  •  SR  (tr  (As PAHR~l  APAfjR~l) 

+  tr(AePAHR-1AuPAHR-1))  (  ' 


Ji„  =  2N  •  R  (tr  (A^PA^R-1  APA^R-1) 

+  tr  (A^PA^R-1  A,,PAhR~1))  *■  ' 


JH  =  2 A  •  R  (tr  (A^PA^R-1  APAf  R-1) 

+  tr  (AePAHR~1AiPAHR~1))  1  ' 


Ju„  =  2N  •  R  (tr  (AwPA"R-[  APA^R-1) 

+  tr(Au,PA"R-lA„PA"R-1))  {  ’ 


Ju(  =  2 N-n  (tr  (AwPAHR-lAPA^R~1) 

+  tr(AuPAHR-lAiPAHR~1))  1  ’ 


=  21V  •  R  (tr  (A^PA^R-1  APAf  R-1) 
+  tr  (A,,PAHR_1  A^PA-^R-1)) 
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